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ABSTRACT 

Transmitted  signals  with  high  time  bandwidth  products  tend  to  resolve 
multiple  reflecting  structural  elements  or  highlights  on  the  body  that  is  being 
illuminated.  This  thesis  develops  a  Kalman  filtering  approach  to  estimating  the 
position  and  velocity  of  the  multiple  highlights  on  a  single  body  undergoing  complex 
motions.  First,  a  general  Kalman  filter  for  direct  recursive  estimation  of  the 
spreading  function  is  derived.  Second,  an  algorithm  which  tracks  peak  locations 
within  the  spreading  function  is  derived  via  an  extended  or  linearized  Kalman  filter. 
The  ability  of  the  Kalman  filter  to  track  kinematic  properties  of  a  multihighlight 
scatterer  is  related  to  the  transmitted  signal’s  mean  squared  bandwidth,  mean  squared 
duration,  and  time-frequency  content  through  the  Cramer-Rao  lower  bound  on 
estimation  errors  for  time  scale  and  time  delay.  It  is  shown  that  the  ability  of  the 
Kalman  filter  to  track  peak  locations  within  the  object  scattering  function  and 
recursively  update  these  peak  locations  depends  strongly  on  the  use  of  signals  with 
high  time-bandwidth  products.  Finally,  a  performance  monitor  which  provides  a 
sound,  monitorable  performance  measure  of  the  tracker  is  introduced  via  the 
innovations  spectrum.  This  performance  monitor  admits  the  ability  of  dynamic 
model  updating  for  adaptive  signal  processing.  This  work  sets  a  groundwork  for 
further  research  into  the  application  areas  of  image  feature  tracking,  robotic  vision, 
high  resolution  radar  and  sonar,  and  medical  imaging. 
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Chapter  1 

INTRODUCTION 


1.1  Motivation 

In  many  applications,  the  propagation  of  a  signal  through  a  medium  is  used  to  derive 
information  about  the  medium  and  about  objects  in  the  medium.  This  method  has  been 
used  in  RADAR,  SONAR,  and  medical  imaging  as  well  a  multitude  of  other  applications. 
A  primary,  often  undesirable  effect  which  manifests  itself  in  this  process  is  the 
phenomenon  of  scattering.  In  the  detection  problem,  (which  is  of  primary  concern  in  this 
thesis)  a  signal  is  transmitted  into  the  medium,  reflects  (scatters)  from  the  object  to  be 
detected  and  is  received  back  at  the  receiver.  The  signal  also  scatters  from  the  other 
objects  in  the  medium,  the  medium  itself,  and  travels  along  multiple  paths  in  the  medium, 
causing  multiple  transmitted  signal  replicas  to  arrive  at  the  receiver  at  different  times 
(delay  spread)[l-3].  For  high  resolution  signals  [4],  the  target  can  be  made  up  of 
multiple  scatterers.  In  this  case  the  signal  may  “over  resolve”  the  target  and  lead  to 
“splitting  loss”  [2].  In  narrowband  systems,  the  problem  of  fading  [2]  is  prominent  for 
distributed  targets  or  spread  channels.  All  of  the  above  effects  give  the  appearance(ffom  a 
detection  point  of  view)  of  multiple  detections.  The  problem  is  then  for  a  signal 
processor  to  determine  which  of  the  returns  (if  any)  are  from  the  target.  The  process  of 
scattering  manifests  itself  in  a  variety  of  forms  including  the  multipath  phenomenon  and 
the  scattering  functions  of  any  objects  in  the  medium  -  including  the  object  to  be 
detected(in  the  case  of  RADAR  and  SONAR). 
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Recently  proposed  [4,6,7]  wideband  processors  such  as  the  wavelet  transform  domain 
(WTD)  estimator-correlator  (EC)  attempt  to  overcome  this  by  modeling  the  multi¬ 
highlight  structure  of  an  object  as  a  spreading  function.  The  multi-path  phenomena  is 
also  treated  as  the  channel  spreading  function.  In  the  narrowband  case,  the  spreading 
function  specifies  how  the  transmitted  signal  is  “spread”  in  delay  and  Doppler.  In  the 
wideband  case,  the  spreading  function  specifies  how  the  transmitted  signal  is  “spread”  in 
delay  and  time-scale.  Spreading  functions  are  two-dimensional  random  functions  and 
their  second  order  statistics  are  given  by  scattering  fimctions[8].  Because  scattering 
functions  can  be  related  to  the  random  scattering  and  propagation  models,  scattering 
functions  are  used  as  prior  statistics  for  implementation  of  model  based  signal  processors 
such  as  the  WTD  EC. 

The  WTD  EC  exploits  the  processing  gain  of  high  resolution  (i.e.  high  time- 
bandwidth  product)  signals  and  uses  an  estimate  of  the  target  and  channel  scattering 
functions  as  prior  statistics  for  implementation  of  the  estimator  branch  of  the  processor 
[9],  The  critical  issue  in  the  performance  of  the  WTD  EC  is  the  accuracy  of  the 
scattering  function  that  is  used  in  the  processor.  A  method  of  updating  the  scattering 
function  estimates  as  new  data  becomes  available  is  required.  To  accomplish  this  update, 
a  Kalman  filter  approach  will  is  used  that  exploits  the  innovations  approach  [10-14]. 

This  method  produces  a  mean-square  sense  optimal  recursive  scattering  function 
estimation  algorithm  and  provides  a  technique  for  monitoring  convergence  of  the 
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estimates.  Furthermore,  the  innovation  method  allows  one  to  establish  the  validity  of  the 
underlying  models. 

1.2  General  Statement  of  the  Problem  -  Objective 

The  main  objective  of  this  thesis  is  to  develop  a  complete  algorithm  for  recursive 
scattering  function  updating  and  estimation.  By  complete,  it  is  meant  that  starting  with 
the  transmitted  signal,  its  shape,  its  parameters  (i.e.  duration  and  bandwidth),  develop  a 
Kalman  filtering  algorithm  which  illustrates  the  effects  of  these  parameters  on  the  filter’s 
ability  to  perform  recursive  estimation.  In  the  course  of  developing  the  Kalman  filter  the 
following  will  be  accomplished: 

1 .  Derive  the  Kalman  filter  for  recursive  scattering  function  updating  and 
estimation. 

2.  Relate  the  Kalman  filter  quantities  to  physical  properties  of  the  signal  and 
scatterer. 

3.  Place  the  Kalman  filter  in  the  context  of  a  model  based  signal  processor,  thus 
illustrating  the  connection  between  tracking  and  detection. 

4.  Using  simulations,  show  the  effects  of  signal  parameters  on  the  ability  of  the 


Kalman  filter  to  track  kinematics. 
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The  first  goal  of  this  thesis  is  to  derive  the  Kalman  filter  for  WTD  scattering  function 
estimates  based  on  the  innovations  approach.  This  derivation  will  show  how  new 
information  (innovation)  will  be  used  to  form  a  new  estimate  of  the  scattering  function 
based  on  the  previous  estimates  and  the  innovations.  It  will  be  necessary  to  derive  a 
WTD  state  equation  for  the  system  under  investigation  as  well  as  a  WTD  measurement 
equation  of  the  spreading  function,  as  this  is  the  quantity  observed.  This  will  be  done 
within  the  framework  of  a  wavelet  transform  based  detection  algorithm. 

The  next  objective  of  the  thesis  is  to  investigate  relationships  between  the  Kalman 
filter  parameters  (e.g.  measurement  error)  and  the  properties  of  the  transmitted  signal. 

The  goal  here  is  to  establish  the  relationship  between  signal  parameters  and  the  Cramer- 
Rao  lower  bound  on  the  error  variance  associated  with  observations  of  time-scale  and 
time-delay. 

The  third  objective  of  the  thesis  is  to  provide  an  implementation  of  the  recursive 
updating  algorithm.  The  results  and  performance  of  the  algorithms  will  then  be  evaluated 
for  selected  examples.  Another  goal  here  is  to  illustrate,  via  the  examples,  the 
relationship  between  the  physical  problem  (target  position  as  a  function  of  time)  and  the 
wavelet  domain  (time/time-scale)  representation  of  the  physical  problem.  Here  the 
method  of  deriving  the  state  equation  and  more  importantly  the  state  transition  matrix 
from  the  physical  problem  is  achieved.  This  critical  step  will  illustrate  the  path  from 
physical  problem  to  scattering  function  prediction  updating  in  the  WTD. 
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Finally,  the  last  objective  is  to  show  simulation  results  which  highlight  the 
transmitted  signal  parameter  effects  and  how  changes  in  the  measurement  error  variance 
(related  to  the  Cramer-Rao  lower  bound)  [15-17]  affect  the  Kalman  Filter.  Here,  a 
possible  performance  monitor  for  the  Kalman  Filter  algorithm,  in  the  form  of  monitoring 
the  “whiteness”  of  the  innovations  sequence[18]  will  be  introduced. 

1.3  Original  Contributions 

The  main  contribution  of  this  thesis  is  the  theoretical  background,  and  algorithm 
development  for  using  a  Kalman  filter  to  track  the  locations  in  the  time  delay  scale  plane 
of  multiple  reflections  which  arise  from  a  single  distributed  abject  which  is  illuminated 
with  a  high  resolution  pulse.  This  differs  from  conventional  tracking/track  association  in 
that  in  these  fields,  each  track  represents  a  different  object,  here,  all  tracks  lie  on  the  same 
object.  This  algorithm  can  be  used  produce  target  spreading  function  estimates  for  use  in 
model  based  signal  processors  such  as  the  wavelet  transform  domain  estimator-correlator. 
In  the  course  of  the  development  of  the  main  theme  of  the  thesis  the  following 
contributions  were  also  made. 

1 .  Derivation  of  the  Cramer-Rao  Lower  Bounds  on  variances  of  estimation  errors  of 
time-delay  and  time-scale  for  a  random  amplitude  scatterer  using  a  wavelet  transform 
detection  scheme. 

2.  Derivation  of  an  extended  Kalman  Filter  for  tracking  kinematic  properties  of  a 
physical  system  using  wideband  probing  signals. 
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3.  Relating  the  parameters  of  the  transmitted  signal  to  the  ability  of  the  Kalman  filter  to 
track  the  kinematics  of  the  physical  system. 

4.  Presentation  of  the  results  of  simulations  which  illustrate  all  of  the  above 
contributions. 

The  main  original  contribution  is  the  development  of  a  unified  approach  to  model- 
based  signal  processing  that  connects  scattering  function  models,  recursive  model 
updating  algorithms  via  Kalman  filtering  theory,  wavelet  transform  implementation  of 
the  EC  and  signal  design  by  CRLB  analysis. 

1.4  Assumptions 

In  this  section,  assumptions  made  in  the  thesis  are  discussed.  The  major  assumptions 
are  the  noises  present  on  the  system  are  uncorrelated  and  Gaussian  and  that  the  scattering 
is  uncorrelated.  This  is  not  always  a  necessary  assumption  and  references  are  provided 
for  those  who  wish  to  extend  the  work  here  to  colored  noises.  Although  multipath 
environments  and  their  properties  are  discussed,  in  the  examples,  the  medium  has  a  single 
propagation  path  that  is  isovelocity.  It  is  assumed  that  the  initial  estimates  are  within  a 
region  of  convergence  for  the  Kalman  filter.  Lastly,  it  is  assumed  assume  in  the 
kinematic  portions  of  the  development  that  the  velocities  are  small  (or  analogously  that 
the  sample  rate  is  high)  so  that  the  overall  data  rate  is  high  enough  for  the  Kalman  filter 
to  track.  This  assumption  is  made  since  the  primary  goal  of  the  thesis  is  to  illustrate  the 
effects  that  the  transmitted  signal  (or  signal  design  issues)  have  on  the  ability  of  the 
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Kalman  filter  to  track.  Also,  spreading  and  scattering  function  models  and  wideband 
representation  (i.e.  time-scale  representation)  of  the  return  are  assumed  to  be  valid.  It  is 
also  assumed  that  all  signals  are  finite  energy,  square  integrable  functions  so  that  they  can 
be  represented  in  the  wavelet  transform  domain.  It  is  also  assumed  that  covariances  of 
received  echos  and  noise  are  positive  definite  functions  that  uniquely  define  a 
reproducing  kernel  Hilbert  space  (RKHS).  Positive  definiteness  is  also  used  to  establish 
existence  of  inverse  kernels  [64],  In  summary,  it  is  assumed  that  all  continuous  signals 
and  stochastic  processes  belong  to  separable  L2  Hilbert  spaces  with  appropriate  inner 
products.  Because  any  infinite  dimensional  separable  Hilbert  space  is  linearly  isometric 
to  I2  continuous  and  discrete  representations  are  equivalent.  Both  representations  are 
used  within  the  thesis  (Section  2.4. 1). 

1.5  Related  Research 

Many  areas  of  application  have  reference  to  tracking,  multiple  target  tracking,  or 
contour  tracking.  In  this  section,  the  goal  is  to  present  related  fundamental  research  in 
the  areas  of  wideband  spreading  and  scattering  and  to  introduce  related  areas  of  current 
research  where  the  features  of  an  “object”  is  investigated.  The  definition  of  “object” 
depends  on  the  application.  First,  the  areas  of  research  which  encompass  wideband 
spreading  and  scattering  are  introduced.  This  research  is  related  to  work  in  the  thesis  in 
that  it  provides  a  starting  theoretical  foundation  on  which  to  derive  the  Kalman  filtering 
algorithms.  Much  of  the  work  done  in  this  area  has  been  to  characterize  time-varying 
propagation  channels  with  multiple  propagation  paths.  This  leads  naturally  a  discussion 
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of  spread  spectrum  systems  and  RAKE  receivers  as  these  systems  attempt  to  recombine 
energy  lost  in  multiple  propagation  environments.  In  this  case  the  “object”  to  be  tracked 
is  the  channel  itself  so  that  an  optimum  recombination  can  be  performed.  The 
relationship  to  the  work  in  the  thesis  is  that  multiple  paths  are  tracked  for  recombination. 
In  multipath  recombination,  the  paths  appear  as  discrete  points  separated  in  delay.  In 
applications  such  as  feature  tracking  in  image  processing,  entire  continuous  curves  which 
appear  in  a  video  stream  are  tracked.  This  has  applications  in  robotic  lip  reading,  or 
tracking  a  person  moving  through  a  crowd:  a  security  application.  Finally,  another  large 
area  of  interest  in  tracking  is  track  association  and  multiple  target  tracking  where  the 
problem  is  to  maintain  tracks  of  multiple  independent  objects.  All  of  these  areas  are 
related  to  the  thesis  through  tracking  theory  (Kalman  filtering).  This  section  points  out 
the  similarities  and  more  importantly  the  differences  between  the  existing  body  of  work 
and  the  work  in  this  thesis. 

In  this  thesis,  the  problem  of  recursively  estimating  the  parameters  of  a  physical 
system  through  the  use  of  probing  signals  is  investigated.  In  proposing  an  algorithm 
using  a  Kalman  filtering  approach,  many  areas  of  current  research  in  the  realm  of 
wideband  signal  processing  [19-24]  are  applicable.  Estimation  as  related  to  narrowband 
spreading  functions  has  received  extensive  coverage  in  the  literature  [16,43],  These 
narrowband  spreading  function  characterizations  have  been  used  for  the  characterization 
of  time-varying  propagation  and  scattering  channels.  Also,  narrowband  source 
localization  and  characterization  has  received  considerable  attention. 
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Characterization  of  these  channels  using  wideband  signals  and  concepts  of  wideband 
spreading  functions  has  received  recent  attention  [24-27],  However,  in  order  to 
understand  the  underlying  principles  of  wideband  spreading,  one  must  also  investigate 
the  properties  of  the  wideband  ambiguity  functions  of  the  transmitted  signal  [28-32], 

The  main  problem  in  this  thesis  is  to  make  recursive  estimates  of  the  properties  of  a 
physical  system  (i.e.,  position  and  velocity)  called  the  state  variables.  One  time 
estimation  of  physical  parameters  has  been  investigated  using  linear  inverse  techniques 
[33],  The  problem  of  identifying  and  localizing  an  object  in  a  changing  environment 
using  wideband  inverse  techniques,  has  also  been  thoroughly  investigated  [24],  The  open 
problem  of  tracking  changes  in  peak  location  within  the  measurement  of  the  wideband 
spreading  function,  and  relating  these  to  the  kinematics  of  an  object  using  a  model  based 
approach  is  investigated  in  this  thesis.  In  the  area  of  model  based  signal  processing  the 
idea  has  been  to  include  models  within  the  signal  processor  to  optimize  system 
performance.  The  work  of  Candy  [88,89]  has  been  to  use  a  model  (environmental)  based 
approach  for  detection  and  localization  of  acoustical  sources.  The  work  in  the  thesis 
differs  in  that  the  models  used  for  tracking  appear  in  the  state  transition  matrix  of  the 
tracker.  These  models  are  of  the  kinematics  of  the  object  which  is  undergoing  rotation, 
translation,  or  both. 


This  formulation  of  the  wideband  spreading  function  estimate  is  based  on  wavelet 
transform  techniques  which  have  received  extensive  recent  attention  [34-40].  In  the  case 
under  investigation,  a  measurement  is  taken  by  using  a  wideband  matched  filter  [4]  which 
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produces  a  range-scale  map  which  indicates  the  “position”  of  the  object  in  time-delay  and 
time-scale.  These  channel  scattering  function  representations  are  exploited  by  spread 
spectrum  communications  for  recombining  the  energy  arriving  at  an  antenna  via  multiple 
propagation  paths. 

In  spread  spectrum  communications  systems,  a  RAKE  receiver  structure  has  been 
used  to  recombine  spread  spectrum  signals  that  propagate  through  slowly  fluctuating 
channels  [76],  In  these  systems,  the  arrivals  of  signals  at  an  antenna  are  only  (in  most 
practical  applications)  separated  in  time  delay  (i.e.,  different  propagation  paths  having 
different  lengths  produce  multiple  arrivals  of  a  transmitted  signal).  The  techniques 
proposed  in  this  thesis  can  be  used  to  extend  this  work  by  making  the  processes 
recursive/adaptive  in  nature  for  two  dimensions.  This  has  application  in  acoustical 
systems  where  time-scale  can  become  appreciable  [76], 

The  closest  areas  of  research  to  which  this  work  can  be  compared  is  in  the  two  areas 
of  contour  and  feature  tracking  in  image  processing  and  in  multiple  target  tracking  and 
track  data  association. 

First  in  the  area  of  image  feature  tracking,  Blake  [77]  provides  an  analog  in  the  sense 
that  two  dimensional  (2D)  image  features  are  tracked  as  they  move  in  an  image  using  a 
model  of  the  motion.  This  is  a  tracking  algorithm  which  predicts  how  features  within  an 
image  will  appear  in  successive  images  based  on  a  model  which  is  used  to  predict  the 
motion  [78-81],  Conditional  density  tracking  uses  learned  dynamical  models,  together 
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with  visual  observations,  to  propagate  a  random  set  over  time.  The  results  of  this  work 
have  been  applied  to  trackers  for  non  rigid  motion  such  as  lip  motion  tracking  for  speech 
recognition.  This  work  differs,  however,  from  the  work  in  the  thesis  in  that  intensity 
values  are  used  to  make  up  the  feature  rather  than  coherent  maps.  Also,  the  work  in  the 
thesis  relates  the  transmitted  signal  parameters  to  measurement  errors  and  ties  these  to 
tracker  performance.  The  measurement  errors  are  less  clearly  defined  in  the  image 
tracking  case  and  vary  from  case  to  case.  Researchers  at  Cameige-Mellon  [82,  83]  have 
used  a  map  correlation  technique  to  track  landmark  features  for  mobile  autonomous  robot 
navigation/vision  systems.  Here  a  normalized  correlation  technique  has  been  used  to 
correlate  masks  of  features  with  masks  extracted  from  robot  images,  the  features  are 
“tracked”  based  on  a  maximized  correlator  output.  Image  feature  tracking,  however, 
does  not  depend  on  transmitted  signals  and  detector  outputs.  Visibility  of  a  feature  in  an 
image  is  analgous  to  detection,  the  difference  lies  in  the  fact  that,  in  this  thesis,  one  of  the 
main  points  is  to  relate  transmitted  signal  parameters  to  tracking  performance.  A  related 
area  of  research  where  the  transmitted  signal  parameters  do  matter  is  in  radar 
tracking/track  association. 

In  the  area  of  radar  tracking,  track  association,  and  data  association  for  multiple  target 
tracking,  much  work  has  been  done  [85-90],  Data/track  association  are  concepts  which 
have  been  used  to  group  measurements  with  the  appropriate  track.  In  the  case  of 
multiple  sensors,  the  data  association  is  done  to  group  the  measurements  from  various 
sensors  with  the  appropriate  target.  Work  has  been  done  to  incorporate  certain  target 
characteristics  into  the  tracking  to  improve  the  data  association  [91].  All  of  this  work. 
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however,  uses  a  separate  tracker  for  each  object  to  be  tracked,  therefore  providing  an 
innovations  sequence  (in  the  case  of  Kalman  filtering)  for  each  target.  In  this  thesis,  a 
single  object,  which  has  multiple  highlights  is  tracked.  The  highlight  motion  is  predicted 
according  to  a  state  transition  which  updates  the  position  of  each  highlight  according  to  a 
kinematic  motion  under  the  assumption  of  rigid  body  motion.  This  method  provides  a 
single  innovations  sequence  which  gives  a  measure  of  tracker  performance.  This  is 
significantly  different  than  using  multiple  trackers  to  track  multiple  objects  and  updating 
each  separately. 

Having  reviewed  the  related  research,  the  next  subject  is  to  determine  the  work  to  be 
done  and  its  relationship  to  the  existing  body  of  related  research.  In  this  thesis,  an 
algorithms  for  tracking  the  peak  values  and  peak  locations  within  a  wideband  spreading 
function  are  derived.  This  extends  the  work  in  the  area  of  spreading  function 
representations  by  giving  a  recursive  spreading  function  estimation  algorithm.  It  will  be 
shown  that  a  model  based  recursive  algorithm  for  tracking  the  kinematic  properties  of 
multiple  highlights  on  a  body  can  be  derived  and  implemented  for  several  models.  This 
extends  work  in  image  feature  tracking  (as  points  on  a  feature  can  be  treated  as  multiple 
highlights).  Also,  work  in  track  association  can  be  extended  by  not  only  tracking  the 
overall  object  trajectory,  but  by  implementing  (for  high  resolution  radar)  a  “sub-tracker” 
which  tracks  the  body’s  highlight  structure.  The  outputs  from  this  filter  can  be  inputs  (in 
the  form  of  range  profiles)  tol  multiple  target  tracking  algorithms  [85,87],  The  following 
section  provides  an  outline  of  the  thesis. 
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1.6  Thesis  Overview 

This  thesis  is  divided  into  six  chapters  including  the  introduction.  The  next  chapter 
contains  an  overview  of  the  theory  necessary  for  completion  of  the  work  in  the  rest  of  the 
thesis.  This  overview  begins  with  a  background  on  Kalman  filtering  and  its  applications. 
The  discussion  then  proceeds  to  cover  the  background  necessary  for  understanding 
wideband  systems,  signals,  and  scattering  phenomena.  Since  a  major  application  of  the 
work  in  the  thesis  is  to  improve  the  performance  of  the  WTD  EC,  the  background 
required  to  understand  and  extend  the  EC  will  be  covered. 

Chapter  3  is  devoted  to  wideband  signal  analysis  in  the  context  of  deriving  the 
Cramer  Rao  lower  bound  on  variances  of  estimation  errors  for  time-delay  and  time-scale 
for  wideband  signals.  Here,  the  derivation  relates  the  lower  bound  on  estimation  error 
variances  to  the  signal  shape  and  signal  parameters.  This  is  important  as  these  quantities 
will  surface  in  the  derivation  of  the  Kalman  filter  in  the  next  chapter. 

Chapter  4  contains  the  area  most  important  in  the  thesis.  Here,  the  WTD  Kalman 
filter  scattering  function  predictor/corrector  is  derived  and  analyzed.  The  parameters  of 
the  wavelet  transform  domain  Kalman  filter  (WTDKF)  are  investigated  in  detail  and 
related  to  the  overall  performance.  These  parameters  are  related  to  those  derived  in 
Chapter  2  (the  CRLB)  and  to  transmitted  signal  parameter. 
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In  chapter  5,  two  examples  are  presented  which  highlight  the  theoretical  results 
derived  in  the  previous  chapters.  The  first  example  is  of  two  spheres  rotating  about  a 
central  point.  The  Kalman  filter  attempts  to  track  the  angle  one  of  the  spheres  makes 
with  the  perpendicular  to  the  line  of  sight,  and  the  normal  velocity  of  the  sphere  as  it 
moves  through  an  arc.  The  second  example  is  a  highlight  moving  in  a  straight  line,  and 
the  filter  attempts  to  track  the  normal  velocity  and  the  time-delay  to  the  scatterer.  This 
chapter  present  results  which  tie  together  the  signal  duration  and  bandwidth,  the  Cramer- 
Rao  lower  bound,  and  the  convergence  of  the  Kalman  filter. 

Chapter  6  contains  the  conclusions  which  may  be  derived  from  the  work  contained  in 
the  thesis.  Here,  suggestions  for  further  research  are  also  presented. 


Chapter  2 

BACKGROUND 


2.1  Introduction 

The  purpose  of  this  chapter  is  to  provide  the  necessary  background  and 
introduction  to  topics  that  are  fundamental  to  the  development  later  in  the  thesis. 

This  chapter  will  provide  an  entry  point  at  which  to  define  notation  used  throughout 
the  thesis.  Finally,  this  chapter  provides  a  discussion  of  previous  research 
fundamental  to  the  development  of  the  work  developed  in  later  chapters. 

Section  2.2  provides  an  introduction  to  the  fundamental  concepts  of  the  wavelet 
transform.  As  the  formulation  of  the  Kalman  filter  is  in  the  wavelet  transform 
domain,  this  background  will  provide  the  necessary  framework  to  understanding  the 
concepts  of  spreading  and  scattering  functions  in  the  wavelet  transform  domain 
(WTD). 

Section  2.3  is  devoted  to  the  concepts  of  the  WTS  wideband  spreading  function 
and  the  second  moment  of  the  spreading  function,  the  real  non-negative  scattering 
function. 


Section  2.4  contains  a  review  of  the  WTD  estimator-correlator  (EC).  This  is 
included  as  the  recursive  scattering  function  updating  algorithm  was  motivated  by  the 
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search  for  an  optimal  recursive  algorithm  for  updating  the  spreading  function  estimate 
in  the  WTD  EC. 

Section  2.5  contains  a  review  of  the  Kalman  filter,  extended  Kalman  filter  (KF) 
and  the  KF  recursive  estimation  algorithm. 

2.1  Wavelet  Transform  Analysis 

This  section  is  to  introduces  the  fundamental  concepts  in  wavelet  analysis  in  the 
context  of  signal  analysis.  It  is  intended  as  an  overview  of  the  concepts  of  time- 
frequency  and  time/time-scale  analysis  and  signal  resolution  properties.  In  this 
section  some  of  the  assumptions  common  to  wavelet  analysis  which  will  impact  later 
discussions  are  reviewed. 

2.1.1  Time-Frequency  Analysis 

The  concept  fundamental  to  the  science  of  signal  processing  for  much  of  the 
last  century  has  been  the  Fourier  transform.  It  has  been  the  basis  for  many  texts  [42]. 
Basically  stated,  the  Fourier  transform  (FT)  is  a  correlation  function.  The  FT 
correlates  a  time  signal  with  a  complex  sinusoid  at  the  desired  frequency  (2.1). 

X(o))  =  F{x(t)}  =  (x(t),e-J6>t\=  ]x{t)e~jmt  dt 

-00 


(2.1) 
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where  <•,•>  represents  an  inner  product  on  the  space  of  square-integrable  functions, 
L2(R).  The  simultaneous  measurements  of  time  and  frequency  are  incompatible  since 
frequency  cannot  be  measured  instantaneously.  Therefore,  (2.1)  represents  the 
average  frequency  content  over  the  length  of  the  signal,  and  gives  us  no  information 
about  when  precisely  the  frequency  occurred  [35], 

To  extend  the  usefulness  of  the  Fourier  transform,  it  can  be  extended  by 
looking  at  a  windowed  portion  of  the  signal  and  applying  a  possibly  complex  window 
g(t)  time-limited  to[-T<t<0]  to  the  signal  and  then  take  the  Fourier  transform. 
Therefore,  the  windowed  signal  contains  only  the  information  in  x(t)  on  the  interval 
[t-T,t].  Taking  the  Fourier  transform  of  the  combined  window-signal  gives  the 
windowed  Fourier  Transform  (WFT). 

xv,  r) = (*„  (o,  )  =  J  msd  -  «■)«-'“'  a  (22) 

—00 

X  contains  both  time  and  frequency  information.  X  contains  information  about  the 
frequency  content  of  x(t)  near  t-x.  The  tradeoff  here  is  that  the  shape,  length,  and 
dimensionality  of  the  window  determines  the  relative  frequency  and  time  resolution 
properties  of  the  WFT.  This  is  known  as  the  uncertainty  principle  [35,42], 
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2.1.2  Time-Scale  Analysis 

Time-Scale  analysis  is  often  associated  with  the  continuous  wavelet  transform 
(CWT).  In  the  CWT  a  time  compression  or  dilation  operation  is  integral  to  the 
process.  The  amount  of  compression  or  dilation  is  referred  to  as  time-scale.  The 
CWT  is  a  generalization  of  the  WFT  in  the  sense  that  the  desired  signal  is  correlated 
with  a  set  of  “wavelets”.  The  CWT  of  a  signal  x(t)  with  respect  to  a  mother  wavelet 
vj/(t)  is  [36] 

X(s,  t)  =  Wvx(s,  r)  =  (x,  y/S  T )  =  J  x(t)  ^*f— 

\  S  -CO  V  S 

One  can  see  that  the  window  in  the  WFT  has  been  replaced  by  the  mother  wavelet  y. 
The  mother  wavelet  must  satisfy  the  admissibility  condition  [36],  Figure  2. 1  is 
provided  to  visualize  the  time  scale  operation. 


Figure  2. 1 :  Morlet  Mother  Wavelet  With  Original  and  Compressed  (s=0.5)  Versions 
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Here  the  relative  delay  has  been  set  to  zero  for  simplicity.  Since  the  compression  of 
the  mother  wavelet  contains  more  cycles  per  unit  time,  it  measures  higher  frequency 
content  in  the  correlation  process.  Dilation  of  the  mother  wavelet  measures  lower 
frequency  content  since  fewer  cycles  are  in  the  same  time  interval. 

Secondly,  an  important  concept  in  wavelet  analysis  is  the  reconstruction  of  the 
original  signal  also  known  as  the  inverse  wavelet  transform  (IWT).  First,  however,  a 
few  conditions  must  be  satisfied. 


If  the  original  signal  is  to  be  reconstructed,  as  said  before  the  mother  wavelet  must 
satisfy  the  admissibility  condition  [36] 


7  |  i//(a>)  |2 

cw  =  2n  \  — ■--■■■  da)  <  oo 

v  L  M 

which  implies  that  the  mother  wavelet  must  have  zero  average  value,  that  it  is  an 
oscillatory  and  bandpass  function. 


(2.4) 


If  the  above  condition  is  satisfied,  the  reconstruction  formula  (IWT)  is  given  by  [36] 

*(0  =  —\\Wy,  x(s,  0 1 ) 

<V  s 


(2.5) 
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2.2.2  Wavelet  Transform  Domain  Spreading/Scattering  Functions 

In  this  section,  the  concepts  of  WTD  spreading  functions  and  the  second 
moment  of  the  spreading  function,  the  scattering  function  are  introduced.  The 
spreading  function  will  be  a  function  of  both  time-scale ,  s,  and  time-delay  ,t.  The 
spreading  function  is  a  set  of  complex  random  variables  which  indicate  how  the 
transmitted  signal  “spreads”  [43,44]  in  scale  and  delay.  The  total  scattering  function 
is  a  group  theoretic  convolution  of  the  scattering  caused  by  individual  scattering 
mechanisms[8, 45,46]. 

The  purpose  of  this  section  is  to  introduce,  qualitatively,  the  physical 
phenomena  of  spreading  functions  and  their  relationship  to  the  quantitative  estimates 
made  on  them  using  the  wavelet  transform. 

In  a  typical  detection  problem,  a  probing  signal  is  transmitted  into  a  medium 
and  reflected  by  motionless  as  well  as  moving  objects  in  the  medium.  These 
reflections  are  received  at  some  receivers)  and  processed.  In  a  complicated 
environment  such  as  the  ocean  and  atmosphere,  the  medium  itself  can  cause  spread  in 
delay  and  Doppler.  This  effect  is  referred  to  as  multipath.  To  illustrate  the 
phenomena,  an  example  of  a  reflecting  object  in  a  multipath  medium  is  given.  If  a 
high  time-bandwith  signal  (see  Chapter  3)  is  used,  the  individual  reflectors  on  the 
target  can  be  resolved.  Therefore,  each  reflecting  feature  causes  an  individual  return 
to  be  “seen”  at  the  receiver.  This  reflector  is  referred  to  as  a  multi-highlight  reflector. 
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The  return  from  the  target  can  be  further  complicated  if  the  medium  contains  a 
multipath.  Figure  2.2  shows  the  geometry  for  a  multi-highlight  object  in  a  multipath 
environment. 


Figure  2.2:  Transmitted  Signal  Scatters  Over  Two  Paths 


Figure  2.2  shows  that  in  this  example,  the  transmitted  signal  reaches  the  target  over 
two  separate  paths.  The  difference  in  length  of  these  two  paths  causes  the  signal  to 
arrive  at  the  target  at  different  times.  Upon  reaching  the  target,  the  signal  is  further 
scattered  by  reflecting  from  each  highlight  on  the  target  as  seen  in  Figure  2.3. 
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Individual  Reflections 


I  xl  |  x2  i  x3  | 

Figure  2.3:  Target  Scatter  for  a  Multi-Highlight  Target 

Now,  mathematical  representations  of  the  above  physical  phenomena  are  given. 

First,  a  model  of  the  echo  received  at  the  receiver  in  Figure  2.2  must  be  introduced. 
This  has  been  discussed  widely  in  the  literature  [26, 33, 47],  The  modeled  echo  is  in 
the  form  of  an  inverse  wavelet  transform 

lit)  =  J  w(x)[U(x)f(t)]dM(x)  +  n(t)  (2«) 

where  w(x)  is  either  the  narrowband  or  wideband  spreading  function,  dp(x)  is  the  left 
invariant  Haar  measure  on  either  the  Heisenberg  (narrowband)  or  Affine  (wideband) 
group  G,  U(x)  is  a  unitary  representation  of  the  corresponding  group,  f(t)  is  the 
transmitted  signal,  and  the  vector  x  represents  either  (x,<J>)  or  (x,j3)  where  t  is  time 
delay,  <|)  is  Doppler  shift,  and  p  is  time  scale.  [U(x)f(t)]  is  given  by  the  following 
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PW(0 ]=r/(»-T)e^'-r> 


BW(0]— £=/1 

■Jfi 


<8 


narrowband 

wideband 


(2.7) 


Now,  in  order  to  obtain  an  estimate  of  the  spreading  function,  the  wavelet  transform 
of  (2. 1)  is  taken  using  the  transmitted  signal  as  the  mother  wavelet,  this  gives 


f  KOPWM]*  dt  =  If  w(x')[U(x)fmU(x')f(t)]’'  dju(x')dt 

.  (2.8) 
+  J  |  n(t)[U(x')f{t)}  dfi(x')dt 


Rewriting  (2.3) 

w(x)=  J \w{x')K{x,x')djLi(x')dt  ^29 

where  K  is  the  auto-ambiguity  function  of  the  transmitted  signal. 

Equation  (2.4)  shows  that  the  estimate  of  the  spreading  function  is  in  the  form  of  a 
convolution  of  the  true  spreading  function  with  the  auto-ambiguity  function  of  the 
transmitted  signal.  If  the  spreading  function  is  within  the  range  space  of  the  wavelet 
transform,  then  K  is  a  reproducing  kernel  for  the  spreading  function.  This  shows  that 
the  properties  of  the  transmitted  signal  have  an  effect  on  estimates  of  the  spreading 


function. 
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An  estimate  of  the  spreading  function  (as  calculated  using  (2.9))  of  the  object 
shown  is  Figure  2.3  is  given  as  the  wavelet  transform  of  the  received  data  under  the 
condition  of  no  noise  and  that  the  reflections  are  mutually  independent.  If  the 
wavelet  transform  of  a  return  from  the  object  in  Figure  2.3  is  taken  four  highlights 
spread  in  delay  and  each  falling  on  a  time-scale  of  1  are  expected.  Figure  2.4  shows 
the  theoretical  wavelet  transform  of  the  return  from  the  object  with  the  transmitted 
signal  as  the  mother  wavelet. 


S=i 


scale 


Delay  (s) 

Figure  2.4:  Object  Spreading  Function  Estimate 

For  this  illustrative  example,  four  copies  of  the  original  signal’s  auto-ambiguity 
function,  spread  in  delay,  and  lined  up  on  a  scale  of  1  are  obtained.  Now,  to  get  the 
total  wideband  spreading  function,  a  group  theoretic  convolution  of  the  target 
spreading  function  with  the  channel  (multipath)  spreading  function  must  be 
performed.  Completing  the  group  theoretic  convolution  for  this  example,  assuming  a 


channel  scatter  as  seen  in  Figure  2.5,  the  overall  spreading  function  estimate  in 
Figure  2.6  is  calculated. 
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Figure  2.5:  Channel  Spreading  Function  Estimate 
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Figure  2.6:  Overall  Spreading  Function  Estimate  for  a  4  Highlight,  2  Multipath, 
Example. 
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As  seen  in  figure  2.6,  fading  due  to  interactions  between  the  highlights  cause  the 
result  that  all  of  the  peaks  are  not  at  the  same  level.  If  there  were  no  fading 
interactions,  eight  distinct  peaks  would  be  present.  Through  the  use  of  high  time- 
bandwidth  (TW)  product  signals,  this  fading  effect  can  be  lessened,  as  the  individual 
peaks  are  sharper  hand  have  smaller  main-lobe  areas  [30]. 

2.2  Estimator-Correlator  Review 

In  this  section,  a  brief  review  of  the  WTD  estimator-correlator  (EC)  [6,7,48,49]  is 
provided.  Although  the  thesis  is  not  about  the  EC,  the  motivation  behind  the  Kalman 
filter  scattering  function  updating  algorithm  lies  within  the  domain  of  the  EC.  A 
major  application  of  the  work  in  the  thesis  is  the  WTD  EC,  and  the  WTD  EC  will  be 
used  in  a  few  of  the  examples  pertaining  to  performance  analysis  of  the  KF. 

In  an  active  system,  a  known  signal,  f(t)  is  transmitted  into  a  channel.  It  is  reflected 
and  received  back  at  the  transducer  as  r(t)  with 

r(t)  =  y(t)  +  n(t)  (2.10) 

Here,  y(t)  is  the  “total  signal”  containing  all  of  the  multipath  and  multihighlight 
components  of  the  echo,  and  n(t)  is  colored  non-stationary  Gaussian  reverberation.  It 
is  assumed  that  the  “total  signal”  y(t)  is  again  given  by  (2.7)[47].  Where  w(x)  is  either 
the  narrowband  or  wideband  spreading  function,  dp(x)  is  the  left  invariant  Haar 


measure  on  either  the  Heisenberg  (narrowband)  or  Affine  (wideband)  group  G,  U(x) 
is  a  unitary  representation  of  the  corresponding  group,  f(t)  is  the  transmitted  signal, 
and  the  vector  x  represents  either  (t,<|))  or  (x,P)  where  x  is  time  delay,  <|>  is  Doppler 
shift,  and  p  is  time  scale.  [U(x)f(t)]  is  again  given  by  the  following 

[£/(x)/(f)]  =  fit  -  t)  eJ^t~x>>  narrowband 

i  (t-A 

[U(x)f(t)]  =  -=f]-j- J  wideband  (2U) 

Equation  (2)  is  a  weighted  integral  of  all  delayed  and  time-scaled  replicas  of  the 
transmitted  signal  with  the  spreading  function  being  the  weight.  This  weighted 
integral  can  be  used  to  model  both  continuously  time-frequency/time-scaled  spread 
signals  or  discretely  distributed  multipaths  and  highlights.  The  total  signal  model 
given  by  (2.6)  can  be  derived  from  the  physics  of  propagation  and  scattering  from 
distributed  targets.  The  delay-Doppler  spreading  function  model  has  been  widely 
used  in  radar  astronomy,  sonar,  and  multipath  communication  [76, 77],  Of  primary 
interest  in  this  work  is  detection  and  characterization  of  scattering  objects  that  have  a 
rich  and  stochastic  highlight  structure.  It  is  also  significant  that  (2.7)  is  in  the  form  of 
an  inverse  time-frequency/time-scale  transform.  This  suggests  that  a  natural 
implementation  of  wideband  EC  is  the  CWTD  implementation.  Next,  the  derivation 
of  the  CWTD  EC  is  reviewed. 

The  EC  computes  the  detection  statistic  /  using 


where  y  is  the  CME  of  y  and  <*,«>r0  denotes  inner  product  on  the  reproducing 
kernel  Hilbert  space  (RKHS)  whose  non-stationary  noise  covariance  function  Ro(t,i) 
is  the  reproducing  kernel  for  that  Hilbert  space. 

In  the  Gaussian  case,  CME  estimates  of  y(t)  can  be  obtained  from  the  orthogonality 
principle  [47] 


<y-y,r>H  =  0 


(2.13) 


where  <»,*>h  denotes  an  inner  product  defined  by 

<x,y>H  =  E[\  x{t)y*  (t)  dt ] 


(2.14) 


For  the  multipath  model  of  (2. 1) 


<y,r>H  =  cf  JJ S(x,x')K(x,x')dn(x')dn(x)  (2.15) 

GG' 


where 


K(x,x')  =  cf{U(x)f,U(x')f)H 


(2,16) 
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is  the  reproducing  kernel  for  the  space  spanned  by  <g,U(x')f>  where  Cf  is  the 
admissibility  constant  defined  by  (2.4). 


Equation  (2. 12)  contains  S(x,x’)  the  inhomogeneous  scattering  function  defined  by 


S(*,x')  =  £{w(*)wV)} 


(2.17) 


The  orthogonality  condition  (2. 10)  itself  is  satisfied  if 


5K0-  /  |5(x,x')(fjK^)0i(v,2)  [t/(*')/(v)]*) 

GG' 

dv  dz[U(x)  f  (t)\dn{x')djj{x) 


where  Qi(v,z)  is  the  inverse  covariance  kernel  for  the  detection  hypothesis  Hi.  That 
is 


J  Qx  (v,  z)Rx  (z,  t)dz  =  S(v  - 1) 


(2.19) 


and 
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Rl(z,t)  =  E[r(z)r\t)]  (2.20) 


2.3.1  EC  -  Continuous  &  Discrete  Representations 

In  this  thesis,  both  the  discrete  and  continuous  representations  of  the  EC  are  used 
equivalently.  In  this  section,  the  relationship  between  the  two  forms  of  the  EC 
detection  statistic  are  reviewed. 


It  has  been  shown  that  from  the  log  likelihood  function  the  detection  statistic  in  the 
form  of  (2.12)  is  obtained  and  that  the  signal  estimate  is  an  integral  equation  (2.6). 
By  vectorizing  the  received  data  r,  signal  s,  and  noise  n  from  equation  (2. 10)  an 
equivalent  discrete  form  of  the  EC  can  be  formed  as  follows: 


Consider  discrete  data. 


r  =  y  +  n  (H^ 

r  =  n  (H0) 


(2.21) 


Where  y  is  the  signal  vector  containing  all  multipath  and  multihighlight  components 
and  n  is  the  noise  vector. 


The  data  covariance  is  given  by 
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E{rrH}  = 


R\  =Ryy  +Rn n 
R0  =  Rnn 


(2.22) 


Assuming  Gaussian  signals  and  noise,  the  detection  statistic  is  proportional  to 
exp(-rHRi'1r)  under  Hi  and  exp(-rHRo'1r)  under  Ho.  Therefore,  the  detection  statistic  is 
given  by 


l  =  rH(Rol-Ri')r 


(2.23) 


This  can  be  written  as 


/  =  r"*f'(*i  -Ro)Ro'r 

=  ((«, -R0)R^r)HR^r 
=  yHRo'r  (2.24) 


Where  y  is  the  conditional  mean  total  signal  estimate.  This  estimate  is  correlated 
with  the  optimally  filtered  data  vector  Ro''r. 
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The  matrix  difference  (Ri-Ro)  is  equivalent  to  the  signal  covariance  Ryy.  Since  the 
signal  and  interference  are  uncorrelated,  Ryy=E{yrH}=RSy  and  Ri=E{rrH}  and  the 
optimal  signal  estimate  can  be  written  as 


R~nr 


(2.25) 


Now  to  get  back  to  the  continuous  representation,  define  the  expansion 


Jt(0  —  xl S  (f) 


(2.26) 


For  a  scalar  process  x(t),  X  is  a  random  vector  and  £  is  a  deterministic  vector  of 
orthonormal  functions  that  make  up  a  complete  orthonormal  set.  The  inverse 
relationship  is  given  by 


x  =  j  x(t)  s  *  (t)  dt 


(2.27) 


Orthonormality  gives 
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\e(u)sH  (w)  du  =  je  (u)s*(u)  du  =  I 

-00  -00 


(2.28) 


Where  I  is  the  identity  matrix.  The  scalar  representation  of  the  signal  vector  is 
therefore 


(2.29) 


By  orthonormality. 


R;!  =  ]e(vy(v)R-'dv 

-00 


where 

00 

r=  J  r{z)e  (z)dz 

-00 


(2.30) 


(2.31) 


Combining  (2.29  ),  (2.30  )  and  (2.31)  gives 


=  ff  (e 1  (0  Ryys'ivMe1  (v)R^e*  ( z))r(z)dxdv 


(2.32) 
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2.3.2  EC  -  WTD  Implementation 

The  WTD  implementation  of  the  EC  follows  if  the  Frobenius-Shur-Godement 
theorem  [36]  is  used  to  compute  the  inner  product  (2. 12).  Equation  (2. 12)  can  also  be 
written  in  the  form  [7] 


l  =  j  Qo  (t,  ■)  r(z)  dzf  dt 


where  Qo  is  the  inverse  covariance  kernel  for  the  detection  hypothesis  Ho  (noise 
only)  and 


ro(0  =  |^o  (Az)  r(z)dz 


(2.34) 


similarly 


>i(0  =  j£?i(C2)  r(z)dz 


(2.35) 


Thus,  the  detection  statistic  can  be  written  as 


(2.36) 
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Then,  according  the  FSG  theorem. 


cf(y,r«)=\{y\u(x)f\)  dfi(x)  (237) 


by  using  equations  (2.34),  (2.35)  and  (2.37),  (2.36)  becomes 


l~Cf\  j S{x,x')  Wfr{x')dn{x')  j  ( Wjr{x))  dju(x)  (2.38) 

G\G' 


where 


wf  r(x)  =  ||  r(z)Qx  (v,  z)[U(x)f(v)]*  dv 


=  \rl(v)[U(x)f(v)Ydv 


(2.39) 


and 


wf  r(x)  =  ||  r(z)Q0(v,  z)[U{x)f{v)]  dv 


=  |^o(v)[C/W/(v)]Vv 


(2.40) 
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For  the  wideband  case,  equations  (2.39)  and  (2.40)  are  the  weighted  wavelet 
transforms  of  r(t)  under  the  appropriate  hypothesis  (corresponding  to  Qo  and  Qi).  The 
final  result  is 


/  =  J  w(x)(pFyV(X))  dju(x) 

G 


(2.41) 


with  admissibility  constant  [36] 


'/ 


-1 


\m\‘ 


dco 


-00 


CD 


(2.42) 


the  mean  square  estimate  of  the  spreading  function  w(x)  given  by 


w(x)  =  J  S(x,x')  W}r(x')d{i(x') 

G 


(2.43) 


This  estimate  of  the  spreading  function  is  important  for  target  characterization. 
Equations  (41)  and  (43)  describe  the  wavelet  transform  domain  implementation  of  the 
EC.  The  processing,  shown  in  Figure  1  is  straight  forward.  The  first  step  is  to 
optimally  filter  [3]  the  data  under  the  Ho  and  Hi  hypothesis.  Then,  using  the 
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transmitted  signal  replica  as  the  mother  wavelet,  correlate  the  wavelet  transform  of 
the  filter  outputs.  This  generates  generalized  range/scale  maps,  of  which  the  Hi 
filtered  map  is  convolved  with  the  generalized  scattering  function  (a  priori  -  from 
models).  The  output  of  this  operation  is  then  correlated  with  the  Ho  filtered  map.  For 
wide-sense  stationary  uncorrelated  scattering  (WSSUS)  the  correlation  is  coherent 
energy  summation.  However,  if  the  WSSUS  assumption  is  not  made,  and  the 
generalized  scattering  function  S(x,x’)  contains  highlight  correlation  information, 
coherent  highlight  integration  may  be  possible. 

2.4  Kalman  Filtering  -  A  Brief  Review 

In  this  section,  the  concepts  involved  in  the  development  of  the  Kalman  filter  are 
investigated.  The  topic  is  approached  from  the  standpoint  of  model  based  signal 
processing.  The  purpose  of  this  section  is  to  provide  the  necessary  background  and 
references  for  the  more  complete  derivation  of  the  recursive  scattering  function 
estimator  developed  later  in  Chapter  4. 

The  problem  that  the  Kalman  filter  has  been  developed  to  solve  is  to  reconstruct 
estimates  of  a  state,  say  x(t),  from  a  set  of  noisy  measurements  y(t)  [51].  There  are 
many  forms  of  the  Kalman  filter  [52-54],  the  predictor-corrector  form  [51]  of  the 
filter  will  be  the  focus.  This  allows  a  model-based  approach  to  the  signal  processing. 
The  major  features  of  the  Kalman  filter  are  presented  in  the  next  section. 
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2.4.1  The  Kalman  Filter  Predictor-Corrector  Algorithm 

The  Kalman  filter  algorithm  requires  the  use  of  several  quantities  that  will  be 
reviewed  throughout  this  section.  Referring  to  the  algorithm  in  Figure  2.7,  individual 
steps  of  the  algorithm  can  be  followed.  First,  at  time  t,  the  previously  filtered 
estimate  x(t-l|t-l)  and  error  covariance  P(t-l|t-l)  exist.  The  goal  is  to  make  the  best 
estimate  of  the  current  state  based  on  the  previous  estimates  (t-1  samples).  This  is  the 
prediction  phase  of  the  algorithm.  Using  the  state-space  model,  the  state  estimate 
x(t|t-l)  and  the  error  covariance  P(t-l|t-l)  are  predicted.  Once  the  prediction  phase  is 
completed,  the  innovations  covariance  R«  and  the  Kalman  gain  K(t)  are  calculated. 
When  a  measurement  y(t)  is  available,  the  innovations,  or  new  information  [1 1,51,55] 
is  determined.  The  old  state  estimate  x(t-l|t-l)  is  used  to  form  the  new,  or  corrected, 
state  x(t|t),  and  covariance  P(t,t).  The  error,  or  innovation,  is  the  difference  between 
the  measurement  y(t|t-l)  and  the  prediction.  The  innovation  is  weighted  by  the 
Kalman  gain  to  correct  the  old  state  estimate  x(t|t-l).  The  error  covariance  is 
similarity  corrected.  Then,  the  algorithm  is  cycled  at  time  t+1.  Figure  2.7  depicts  the 
entire  algorithm  [51]. 


39 


Figure  2.7  Kalman  Filtering  Algorithm 

Candy  [51]  provides  an  example  of  the  algorithm  through  two  time  steps  which 
provides  insight  into  the  operation  of  the  filter  and  the  calculation  of  the  individual 
quantities.  The  example  is  reproduced  here.. 
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Example  4. 1  This  example  illustrates  the  calculations  performed  in  the  algorithm  in 
the  above  section  given  initial  estimates.  The  calculations  are  carried  out  for  times 
t=l  and  t=2.  State  noise  in  this  example  is  zero. 


t  =  1 


Start  Time 


*(1|0)  =  ^(0)*(0,0) 

£(1|0)  =  ^(0)£(0|0K(0) 

Re  (1) =  C(l)^(l  1 0)C'(1)  +  Rv  (1) 

m  =  P0.\0)C'(\)R;l(\) 
e(l)  =  y(l)-C(l)x(\\0) 
jc(l|l)  =  Jc(l|0)  +  A’(l_e(l) 

P(1|1)  =  [/-^(1)C(1)]P(1|0) 

t  =  2 

x(2 1 1)  =  ^(l)jc(l,l) 

P(2|1)  =  ^(1)P(1|1M'(1) 

Re  (2)  =  C(2)P(2 1 1)C'(2)  +  Rv  (2) 
K(2)  =  P(2\\)C’(2)R;\2) 
e(2)  =  y(2)-C(2)x(2\\) 
x(2 1  2)  =  jc(2  1 1)  +  K{2)e{2) 
P(2\2)  =  [I  -  K{2)C(2)]P{2 1 1) 


Estimate  x(l|0)  using  state  transition  A 

Estimate  error  covariance  P  using  state 
transition  A 

Calculate  Innovations  Covariance  at  t=l 
Using  Measurement  Matrix  C 

Calculate  Kalman  Gain  K  at  t=l 

Use  Measurement  y  at  t=l  to  calculate 
innovations  e 

Calculate  Updated  State  Estimate  at  t=2 
Calculate  Error  Covariance  Estimate  at  t=2 


Increment  Time 


Estimate  x(2|l)  using  state  transition  A 

Estimate  error  covariance  P  using  state 
transition  A 

Calculate  Innovations  Covariance  at  t=2 
Using  Measurement  Matrix  C 

Calculate  Kalman  Gain  K  at  t=2 

Use  Measurement  y  at  t=2  to  calculate 
innovations  e 

Calculate  Updated  State  Estimate  at  t=2 
Calculate  Error  Covariance  Estimate  at  t=2 


This  is  simply  a  brief  review  of  the  operation  of  the  Kalman  filter.  Chapter  4 
provides  a  detailed  derivation  of  the  Kalman  filtering  algorithm  used  as  well  as  a 
detailed  analysis  of  all  of  the  quantities  involved.  In  this  problem  tracking  of  the 
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Kinematic  properties  of  a  system  using  probing  signals  as  the  measurements  of  the 
state  will  be  done.  The  state  of  the  system  will  be  derived  as  well  the  measurement 
equation.  The  main  focus  will  be  on  the  problem  of  designing  a  signal  that  provides 
the  resolution  necessary  for  the  Kalman  filter  to  track.  Also  discussed  is  how 
quantities  which  arise  within  the  Kalman  filter  derivation  can  be  used  to  provide 
performance  monitoring  of  the  entire  system. 

2.5  Application  of  the  Target  Scattering  Function  Algorithm  in  the  WTD  EC 

In  this  section  the  applicability  of  the  Kalman  filtering  algorithm  in  an  adaptive 
model  based  signal  processor,  namely,  the  WTD  EC  is  discussed.  The  WTD  EC  uses 
an  estimate  of  the  wideband  spreading  function  of  the  target  (and  channel)  in  the 
calculation  of  the  EC  detection  statistic.  The  output  of  our  Kalman  filter  will  be  an 
estimate  of  the  target  portion  of  the  spreading  function  for  use  in  the  estimator  branch 
of  the  EC.  This  is  illustrated  in  Figure  2.8. 
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Figure  2.8:  Application  of  Kalman  Filter  to  WTD  EC 


In  this  case,  an  attempt  to  update  the  target  spreading  function  estimate  based  on  the 
kinematics  of  the  problem  would  be  done.  The  EC  detector  output  would  be 
maximized  on  a  ping  to  ping  basis,  and  provide  performance  monitoring  of  the 
overall  system. 


Chapter  3 


ESTIMATION  ERRORS  -  THE  CRAMER-RAO  LOWER  BOUND 

3.1  Introduction 

In  this  chapter,  the  Cramer-Rao  (C-R)  bound  [16,17]  for  the  minimum  variance  of 

joint  estimates  on  time-scale  and  time-delay  is  derived  for  wideband  echoes  that  are 
processed  with  a  wideband  (wavelet  transform)  matched  filter.  The  C-R  bound  on 
Doppler  shift  for  narrowband  signals  is  well  known  [16],  but  the  extension  to  the  more 
general  case  of  high  time-bandwidth  signals  and  time-scale  is  not.  The  typical  measure 
of  resolution  of  a  signal  is  the  width  at  the  -3dB  points  on  a  transmitted  signal’s  auto¬ 
ambiguity  function.  This  measure  is  meaningful,  however,  only  in  the  absence  of  noise. 
When  noise  of  any  kind  is  present,  the  resolution  of  the  signal  is  limited  not  only  by  the 
main  lobe  width  (a  function  of  time  bandwidth  product  (TW)),  but  also  by  the  signal  to 
noise  ratio.  Therefore,  in  order  to  accurately  speak  of  resolution,  the  C-R  bound  must  be 
used,  as  it  is  a  measure  of  the  coupled  effect  of  noise  and  main  lobe  width. 

3.2  The  Auto-Ambiguity  Function  -  Resolution 

Resolution  is  the  property  by  which  one  can  distinguish  (in  both  time-scale  and  time- 
delay)  between  separate  returns  in  a  range-scale  map.  Since  in  this  thesis,  the  detection 
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and  estimation  problem  is  under  consideration,  this  is  a  good  place  to  lay  out  the  problem 
in  context. 

3.3  The  Detection  and  Estimation  Problem  -  No  Noise 

In  this  problem,  a  high  resolution  signal  (large  TW)  is  transmitted  into  a  medium  in 
which  it  is  desired  to  resolve  two  independent  point  scatterers  which  are  separated  by  a 
distance  d  (Figure  3.1)  moving  with  normal  velocities  vl  and  v2. 


Figure  3.1:  Problem  Geometry 

As  seen  in  Figure  3. 1,  the  first  object  is  a  distance  r  from  the  transmitter/receiver  (Tx/Rx). 
Both  objects  are  moving  with  normal  (toward  or  away  from  the  Tx/Rx)  velocities  vl  and 
v2.  A  high  TW  signal  is  transmitted  into  the  medium  (homogenous/infinite)  and  reflects 
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from  the  objects.  This  causes  the  total  echo  to  contain  two  echoes  each  of  which  is  time- 
scaled  (compressed  or  expanded  depending  on  the  direction  of  the  motion).  Now,  the 
data  is  match-filtered  using  the  transmitted  signal  as  the  mother  wavelet,  a  range  -  scale 
map  of  the  echo  (Figure  3.2)  is  produced. 


Time-Scale 


Range (m) 


Figure  3.2:  Range-Scale  Map 


where  each  of  the  peaks  in  the  figure  represents  the  delay  and  time-scale  location  of  the 
individual  scatterers.  Now  if  the  signal  has  less  resolution  and  the  scatterers  are  close 
together,  as  is  seen  in  Figure  3.3,  the  location  of  the  actual  peak  becomes  increasingly 
ambiguous  as  the  signal  ambiguity  function  becomes  less  sharp. 


Delay 

Figure  3.3:  Ambiguity  Main  Lobes  Begin  To  Overlap 

As  seen  in  Figure  3.3,  when  the  resolution  is  low,  the  ability  to  distinguish  individual 
scatterers  in  the  range-scale  plane  is  reduced.  In  Figure  3.4,  it  is  seen  that  two  scatterers 
that  are  closer  in  range  are  indistinguishable  as  the  ambiguity  functions  overlap. 


Figure  3.4:  Scatterers  Become  Indistinguishable 
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In  summary,  the  higher  the  TW  product,  the  smaller  the  main  lobe  area  [6],  thus 
increasing  the  ability  to  distinguish  individual  scatterers.  The  addition  of  noise  in  any 
form  further  complicates  the  ability  to  estimate  the  position  of  the  scatterers  in  the  time¬ 
scale  plane. 

3.4  Detection  and  Estimation  -  Additive  Noise 

Now  the  effects  of  adding  noise  of  any  kind  to  the  problem  are  investigated.  Plotting 
the  ambiguity  function  for  multiple  delay  hypothesis  and  for  a  single  scale  hypothesis  so, 
a  plot  seen  in  Figure  3.5  is  obtained. 


Scale 


Figure  3.5:  Ambiguity  Function  Without  Noise  for  Scale  Hypothesis  so 
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Adding  noise  to  the  time-series  that  produced  the  ambiguity  function  has  the  effect  that 
the  actual  delay  (i.e.,  the  delay  to)  is  much  more  difficult  to  estimate  (Figure  3.6).  The 
estimated  delay  in  not  equal  to  the  actual  peak.  This  effect  is  exaggerated  at  lower 
signal-to-noise  ratios  and  broader  ambiguity  function  main  lobes.  It  is  the  bound  on  the 
variance  of  this  estimate  that  is  the  desired  quantity.  This  bound  is  called  the  Cramer- 
Rao  lower  bound  [56], 


x  estimated 


Scale 


Figure  3.6: 

In  summary,  it  is  desired  to  calculate  a  bound  on  the  variance  of  estimates  that  can  be 
made  on  the  location  of  a  scatterer  in  the  scale-delay  plane.  This  bound,  the  C-R  bound. 


To  Delay 

Ambiguity  Function  With  Noise  Added 
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will  include  the  coupled  effects  of  ambiguity  main  lobe  area  and  signal-to-noise  ratio. 
This  is  a  useful  quantity  for  both  signal  design  and  detector  performance  analysis. 


3.5  C-R  Lower  Bound  Derivation 

In  this  section  the  C-R  bound  starting  with  the  general  wideband  representation  of  a 
received  echo  with  noise  added  is  derived.  In  order  to  simplify  the  derivation  it  is 
assumed  that  the  noise  is  uncorrelated  with  the  signal,  that  it  is  white,  and  is  Gaussian 
distributed. 

3.5.1  Wideband  Representation  of  the  Echo 

Assuming  that  the  scatterer  is  a  slowly  fluctuating  point  target  moving  with 
constant  velocity  through  a  homogenous,  infinite  medium.  The  echo  from  the 
scatterer  received  at  the  receiver  (Figure  3. 1)  is  modeled  as  [4,6,7,16] 


Ht)  =  JF,b -fs/(s(t  -  r»  +  n(t) 


(3.1) 


here  Et  is  the  transmitted  signal  energy,  b  is  a  zero  mean,  complex,  Gaussian 
random  variable,  and  E{|b2|}=2cb2.  n(t)  is  an  independent,  zero-mean,  white, 
Gaussian  process  with  covariance  given  by 
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E[n(t)n(t)]  =  N08(t-u)  (32) 

From  Equation  (3. 1),  s  and  x  are  the  non-random  scale  and  delay  parameters  to  be 
estimated  and  for  which  the  C-R  bound  will  be  calculated. 

3.5.2  The  Detector 

In  order  to  continue  the  derivation,  it  is  necessary  to  define  the  detector  structure 
being  used  (i.e.,  define  the  detection  statistic  as  a  function  of  s  and  x).  In  this 
case,  a  wideband  matched  filter  operation  (wavelet  transform)  given  by 

L(s,  t)  =  |  r(tyjsf*  (s{t  -  r))dt 

—00 

(3.3) 


is  performed.  Where  L  is  the  detection  statistic,  s  is  the  hypothesized  time  scale,  x 
is  the  hypothesized  time  delay  and  r(t)  is  given  by  (3.1).  The  next  step  is  the 
derivation  is  to  derive  the  likelihood  function.  This  has  been  done  in  [16]  and  the 
result  stated  here 


In  A(s,  r) 


1  E, 
N0N0+Er 


(3.4) 


Here  No  is  the  noise  energy  and  Er  is  the  total  echo  energy. 
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Having  obtained  equations  (3.3)  and  (3.4),  the  next  step  is  to  investigate  the 
Fisher  information  matrix  [57,58], 

3.5.3  The  Fisher  Information  Matrix 

In  order  to  estimate  the  non-random  parameters  s  and  t,  the  Fisher  information 
matrix  must  first  be  derived.  The  lower  bounds  on  the  estimates  (C-R  bounds)  are 
given  by  the  diagonal  components  of  the  inverse  of  the  Fisher  information  matrix. 
The  Fisher  information  matrix  is  derived  in  detail  in  [16]  the  results  are  stated 
here.  For  a  given  likelihood  function  the  elements  of  the  information  matrix  are 
given  by  the  following 


d2  In  A 

ddjdcij 


(3.5) 


Where  InA  is  given  by  (3.4)  and  a;  and  aj  are  s  and  t  respectively.  As  these  are 
non-random  parameters,  the  expectation  in  (3.5)  is  over  r(t)  and  n(t).  Using  (3.5) 
the  individual  elements  of  the  information  matrix  are  given  by  [56,57] 


d2  InA 
dr1 


J\\  -~E 


(3.6) 
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d2  In  A 
~d?~ 


(3.7) 


J 12  ™"  J  w  ~  ~~E 


21 


d2  InA 
dsdr 


(3.8) 


Inserting  (3.4)  into  (3.6), (3. 7),  and  (3.8)  produces 


Jn=-E 


2c  Re 


dLdL*  _  d2L* 

- +  l — 

dr  dr  dr2 


(3.9) 


J  22  =  ~E 


2cRe 


dLdL 
ds  ds 


+  L 


d2L* 

ds2 


(3.10) 


2c  Re 


az,  azT 

ds  dr 


+  L 


dsdr 


(3.11) 


where 
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C  = - L-=^ 

N o  N0  +  E0 

It  is  now  necessary  to  evaluate  each  of  the  equations  (3.9),  (3. 10),  and  (3.11) 
separately.  As  the  algebra  is  intense,  only  the  term  Jn  will  be  derived  (Appendix 
A)  in  detail.  The  other  terms  have  very  similar  derivations. 

The  individual  terms  as  given  by  (3.9),  (3. 10),  and  (3.1 1)  with  (3.3)  substituted 
for  L  are  given  in  final  form  as: 


Jn  =  ^ 


f  6)'' 


(3.13) 


J  22  —  2  cEr 


4  ft2Ta>t)-  —  +  at 
v  ’  4  r 


(3.14) 


As  can  be  observed  in  equation  (3.13),  the  Jn  term  is  a  function  of  mean  squared 
bandwidth  indicated  by  the  bar  and  the  hi  term  is  proportional  to  a  quantity  a, , 
called  the  alpha  moment  [30],  By  inverting  these  quantities,  the  C-R  bounds  on 
time-delay  and  time-scale  are  calculated.  In  the  next  section,  illustrative 
examples  of  these  quantities  as  a  function  of  signal  duration,  bandwidth,  and 
signal  to  noise  ratio  are  given. 
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3.6  Examples  of  the  C-R  Lower  Bound 

In  this  section,  example  calculations  of  the  CRLB  for  time  delay  and  time  scale  are 
presented.  In  order  to  see  the  effects  of  changing  various  transmitted  signal  parameters,  a 
few  simplifying  assumptions  are  made  here.  The  first  major  assumption  is  that  the 
estimation  errors  of  time-scale  and  time-delay  are  uncoupled  [16].  In  order  for  this  to  be 
true,  the  signal  envelope  must  be  real  and  the  signal  must  be  symmetric  about  the 
frequency  axis  [59],  Making  these  assumptions,  only  the  Jn  and  J22  terms  are  needed  to 
calculate  the  C-R  bounds.  This  simplifies  the  calculations  as  the  cross  terms  J12  are  not 
needed  any  longer.  It  is  a  research  issue  to  investigate  further  the  signal  design  problem 
in  terms  of  estimation  error  decoupling.  The  C-R  bounds,  making  the  above 
assumptions,  can  be  written  as  [60]: 


var[r  -  f]  =  (2  cEr ) 


-1 


co 

SJ 


V 


00 

Vsj 


-1 


(3.15) 


And 


var[s-s]  =  (2cEr) 


-1 


4x2  \coi)-^+at 


-1 


(3.16) 


As  expected,  the  bounds  are  functions  of  signal  parameters  as  well  as  signal  to  noise 
ratio.  The  following  sections  contain  examples  of  the  CRLB  for  time-scale  and  time- 
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delay  to  illustrate  the  dependance  on  signal  parameters.  Also,  it  is  well  known  that  the 
alpha  moment  appearing  in  equation  (3.16)  above  is  not  only  a  function  of  signal 
parameters,  but  also  a  function  of  the  signal  shape.  The  alpha  moment  is  larger  for 
signals  which  have  their  higher  instantaneous  frequency  content  located  at  the  end  points 
of  the  signal  [61,62],  This  has  an  effect  on  the  CRLB,  and  therefore,  two  signal  types 
will  be  used  to  illustrate  the  point.  A  VFM  signal  which  has  a  time-frequency  plot  as 
given  in  Figure  3.7  and  a  roof-top  FM  (RTFM)  which  has  a  time-frequency  plot  as 
plotted  in  Figure  3.8  will  be  used  for  the  examples  to  follow. 
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Figure  3.7:  Time  Frequency  Plot  for  VFM  Signal 


56 


2000 

1950 

1900 

1850 

Frequency  l8W 
(Hz)  1750 

1700 

1050 

1600 

1550 

1500 

-0.5  -0.4  -0.3  -0.2  -0.1  0  0.1  0.2  0.3  0.4  0.5 

Time 

Figure  3.8:  Time  Frequency  Plot  for  Roof-Top  Signal 

3.6.1  CRLB  for  Time  Delay 

In  this  section  plots  of  the  CRLB  for  time  delay  are  given.  These  plots  are  only  given 
for  the  VFM  signal  as  the  CRLB  for  time  delay  is  independent  of  signal  shape.  It  is  also 
seen  that  the  CRLB  for  time-delay  is  a  function  only  of  mean  squared  bandwidth  and 
signal  to  noise  ratio.  Figure  3.8  is  a  plot  of  the  CRLB  for  time  delay  as  a  function  of 
mean  squared  bandwidth. 
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Figure  3.9:  CRLB  for  Time  Delay  for  Varying  SNR  as  a  Function  of  Bandwidth 

Figure  3.9  confirms  the  expectation  raised  in  section  3.2,  that  is,  as  the  SNR  is  decreased, 
progressively  worse  estimates  on  the  time  delay  are  made.  Also,  increased  bandwidth 
can  recover  some  of  the  loss  due  to  signal  to  noise  ratio,  this  is  a  positive  argument  for 
using  wideband  signals  for  estimation  of  range.  Figure  3.9  shows  that  as  the  bandwidth 
is  increased,  the  CRLB  is  lowered.  This  trend  shows  that  for  a  fixed  signal  to  noise  ratio, 
statistically  speaking,  by  using  signals  with  the  higher  frequency  components  at  the  ends 
of  the  signal  in  time,  better  estimates  of  range  can  be  made  in  the  detection  problem. 

How  these  progressively  worse  estimates  affect  the  ability  of  the  Kalman  filtering  routine 
to  converge  is  the  subject  of  Chapter  5. 
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3.6.2  CRLB  for  Time-Scale 

In  this  section  the  more  complex  issue  of  the  CRLB  for  time-scale  estimation  errors  is 
investigated.  In  this  case,  the  CRLB  is  not  only  determined  by  the  transmitted  signal 
parameters  such  as  duration,  bandwidth,  and  center  frequency,  but  also  on  the  alpha 
moment,  which  is  itself  a  function  of  the  signal  shape  [61-63],  This  section  will  present 
examples  that  illustrate  this  fact.  In  this  section,  the  CRLB  for  time-scale  will  be 
calculated  for  a  variety  of  transmitted  signal  parameters  and  for  the  two  signals 
envelopes:  VFM  and  RTFM.  It  is  a  simple  matter  to  derive  the  alpha  moment 
analytically  for  the  VFM  and  RTFM  using  the  equation  [60, 61, 62] 

T_ 

at=—  \t2co2{f)dt  (3- 17) 

T  _r 

~2 

If  the  inverse  of  the  alpha  moment  (which  is  proportional  to  CRLB)  it  is  seen  that,  as 
expected,  the  CRLB  for  time  scale  is  lower  when  the  instantaneous  frequency  content  is 
higher  at  the  endpoints  of  the  signal.  Figure  3. 10  is  a  plot  of  the  CRLB  for  time-scale  as 
a  function  of  duration  for  the  VFM  and  RTFM  for  equal  mean  square  bandwidths  of  1000 


Hz. 
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Figure  3. 10:  CRLB  for  VFM  and  RTFM  vs.  Duration 

As  can  be  seen  in  Figure  3.10,  the  CRLB  is  lower  for  the  VFM  for  the  same  mean 
squared  bandwidth  and  durations.  This  indicates  that  better  estimates  can  be  made  if  the 
higher  frequency  components  are  pushed  out  to  the  ends  of  the  signal  (in  time). 
However,  as  is  seen  in  Figure  3.1 1,  this  effect  is  lessened  as  the  total  bandwidth  gets 
much  higher. 
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Figure  3.11:  Time  Scale  CRLB  for  5kHz  Bandwidth  VFM  and  RTFM 

3.7  Conclusions 

The  CRLB  for  time-delay  and  time-scale  for  the  case  of  wideband  signals  has  been 
derived  and  their  dependence  on  the  transmitted  signal’s  parameters  have  been 
investigated.  As  expected,  it  has  been  shown  that  as  the  effective  duration  is  increased, 
the  ability  to  estimate  time-scale  is  increased.  Also,  as  the  mean  squared  bandwidth  is 
increased,  the  ability  to  estimate  time-delay  is  improved.  It  was  shown  that  the  CRLB  for 
delay  and  scale  can  be  related  to  the  parameters  of  the  transmitted  signal  and  to  the  alpha 
moment.  The  alpha  moment  is  a  function  of  signal  bandwidth  as  well  as  signal  shape. 
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The  important  conclusion  is  that  the  CRLB  for  the  variance  on  estimation  errors  made  on 
time-scale  and  time-delay  can  be  related  to  the  time-bandwidth  product  of  the  transmitted 
signal.  This  is  an  important  signal  processing  issue  as  it  has  been  shown  that  processing 
gain  (in  the  detection  problem)  is  proportional  to  time-bandwidth  product  for  point 
scatterers.  The  purpose  of  this  derivation  and  analysis  in  terms  of  the  overall  thesis  is 
that  the  CRLB  for  delay  and  scale  dictates  a  minimum  variance  on  the  errors  associated 
with  estimating  time-delay  and  time-scale.  These  quantities  will  be  used  in  the  next 
chapter  in  the  form  of  a  measurement  equation  in  the  Kalman  Filter.  In  the  forthcoming 
development,  the  performance  of  the  Kalman  filter  will  be  shown  to  be  tightly  tied  to  the 
measurement  error  covariance  matrix.  As  it  will  be  shown,  the  elements  within  the 
measurement  error  covariance  matrix  are  the  variances  predicted  by  the  CRLB  .  This 
relationship  is  critical,  as  the  best  possible  performance  achievable  in  the  Kalman  filter 
depends  on  (among  other  things)  the  CRLB  for  time-scale  and  time-delay.  Having 
shown  that  the  CRLB  is  tied  to  the  transmitted  signal  parameters,  a  complete  concept  can 
be  formed,  tying  realistic  measurement  quantities  to  expected  Kalman  filter  estimator 
performance. 


Chapter  4 


Wavelet  Domain  Scattering  Function  Predictor/Corrector  Kalman 

Filtering 


4.1  Introduction 

As  discussed  in  Chapter  1,  two  approaches  to  updating  the  a  priori  scattering 
function  information  are  investigated  in  the  thesis.  The  first,  to  be  discussed  in  this 
chapter  is  the  direct  recursive  estimation  of  the  spreading  function  based  on  spreading 
function  estimates  using  a  Kalman  filter.  Upon  completion  of  the  derivation  it  was 
found  that  the  complexity  of  the  results  (i.e.,  four  dimensional  state  transition  and 
measurement  equations)  pointed  to  approaching  the  problem  using  a  simplified 
approach.  The  simplified  approach  is  the  tracking  of  peak  locations  within  the 
spreading  function  using  an  extended  Kalman  filter.  This  approach  is  more  easily 
implemented,  and  can  be  used  to  update  spreading  function  estimates  by  updating  the 
scattering  function  found  in  Equation  (2.40).  The  simplified  approach  is  introduced 
at  the  end  of  Chapter  4,  and  in  implemented  in  Chapter  5. 

This  chapter  contains  the  derivation  and  analysis  of  the  Kalman  filtering  approach 
to  scattering  function  recursive  estimation  and  updating.  Presented  here  are  the 
results  of  the  derivation  of  the  Kalman  filter  for  WTD  scattering  function  estimation 
and  updating  based  on  the  innovations  approach  [1 1,5 1].  The  parameters  of  the 
Kalman  filter  (e.g.,  the  criterion  and  quality)  are  analyzed  in  terms  of  the  WTD 
observables  (e.g.,  spreading  and  scattering  functions)  and  the  performance  parameters 
of  the  Kalman  filter  are  compared  to  and  analyzed  with  the  performance  measures  of 
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the  WTD  EC.  It  is  shown  that  relationships  exist  between  the  scalar  figure  of  merit 
derived  in  Chapter  2  and  the  criterion  of  the  Kalman  filter  (i.e.  the  trace  of  the  error 
covariance).  It  is  also  shown  that  the  EC  detection  statistic  error  function  is  related  to 
the  Kalman  filter  error  covariance. 

4.2  Innovations  Approach  to  the  Kalman  Filter 

In  this  section,  the  Kalman  filter  for  recursive  scattering  function  updating  using  the 
innovations  approach  outlined  by  Kailath  [1 1]  is  derived.  The  general  introductory 
theory  of  the  Kalman  filter  operation  was  given  in  Section  2.2. 

4.2. 1  State  Equation 

The  first  relationship  required  in  the  derivation  of  the  Kalman  filter  is  referred 
to  as  the  state  equation  [51].  Here,  a  state  equation  in  terms  of  the  spreading 
function  will  be  given  as 


w(x)  =  A(x  -  1)w(jc  - 1)  +  u(x  - 1)  (4.1) 

Here,  x  is  given  as  X€R2  or  x=[t  s]  where  x  is  the  variable  of  time  delay  and  s 
is  the  variable  of  time-scale.  The  state  equation  gives  the  relationship 
between  the  x-1  (previous  scale  and  delay)  and  the  current  scale  and  delay  x. 
The  relationship  (1)  therefore,  gives  the  relationship  of  the  new  spreading 
function  w(x)  to  the  previous  spreading  function  w(x-l)  where  the  (x-1) 
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portion  of  the  expression  indicates  the  previous  s  and  t  (i.e.  if  x=[s  t]  then  x- 
l=[st-i  tt-i].  This  is  accomplished  by  a  known  variable  A  called  the  state 
transition  matrix.  The  state  transition  matrix  is  an  appropriately  dimensioned 
matrix  which  is  derived  from  the  physics  of  the  problem  (see  Section  5.2  and 
5.3  for  specific  examples).  For  the  purposes  of  the  derivation  it  is  sufficient  to 
state  that  the  state  equation  (1)  models  the  system  under  investigation  in 
general.  A  is  an  appropriately  dimensioned  matrix  called  the  state  transition 
matrix.  Here,  q  is  the  state  noise  which  is  assumed  to  be  zero  mean  and  white 
with  covariance  Rq  and  q  is  uncorrelated  with  w,  the  spreading  function. 


The  next  step  in  the  derivation  is  to  arrive  at  a  measurement  equation. 


4.2.2  The  Measurement  Equation 

In  order  to  begin  the  Kalman  filter  derivation,  it  is  first  necessary  to  find  (for  our 
EC  processing  structure)  the  proper  representation  of  the  measurement  equation  of 
the  general  form 


>%(*)  =  C(*K(*)  +  «(*)  (4.2) 

here  wa  is  the  observed  spreading  function,  C  is  a  matrix  operating  on  the  desired 
variable,  and  n  is  a  noise  process.  The  subscript  m  denotes  the  fact  that  this  is  an 
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observed  (measured)  variable  and  the  subscript  o  denotes  that  this  is  the  desired 
variable  or  the  variable  to  be  estimated  from  the  observations. 


As  seen  in  Section  2.3  in  the  outline  of  the  WTD  EC,  the  measurement 
equation  under  reverberation  limited  conditions  is  given  by: 


r(t)  =  y(t)  +  n(t)  (4  3 

where  r(t)  is  the  received  data,  y(t)is  the  signal  with  reverberation ,  and  n  is  a  colored 
noise  process.  Recall  Section  2. 1  that  y(t)  for  the  EC  is  given  by: 


m = J  <*)[c/(*)/(f)]rfM*) 

G 


(4.4) 


where  w(x)  is  the  spreading  function  of  Section  2.3,  U(x)  is  a  unitary  operator 
operating  on  either  the  Heisenberg  (narrowband)  or  affine  (wideband)  group,  f(t)  is 
the  transmitted  signal  and  p(x)  is  the  left  invariant  Haar  measure  for  the  appropriate 
group. 

Substituting  (4.4)  into  (4.3)  and  taking  the  wavelet  transform  of  both  sides: 

J  r(0[tw«)]’<*  =  J  (  »'(*')  [tf  (*)/(/))*  du(x')  dt  +  f  n(t)[U(x)f(t)]’dt  (4  5) 


In  the  first  integral  on  the  right  hand  side  the  time  dependent  terms  can  be  grouped 
together  in  the  time  integral  pulled  into  the  G  integral  resulting  in 
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\ r(l)[U(x)f =  J  w(x')  ([[£/(*)/«]  [t/W/«r  dt  J i^x)  +  J «(/)[£/«/(, )]V< <4-6) 


Recognizing  that  the  time  integration  in  the  first  term  on  the  right  hand  side  of  (4.6) 
represents  the  reproducing  Kernel  which  spans  the  space  of  signals  admissible  f(t), 
(4.6)  can  be  written  as 


"mi*)  =  jw0(x')  K(x,x')  dM(x')  +  \n{t){U{x)f{t)Ut 


Noting  that  the  second  term  is  the  wavelet  transform  of  the  received  noise,  for  brevity 


=  d/u(x)  +  n^(x) 


where  the  subscript  wf  denotes  the  wavelet  transform  with  f(t)  as  the  mother  wavelet. 
The  measurement  equation  is  beginning  to  take  shape.  It  is  now  necessary  to  show 
that  the  integral  operator  (really  a  convolution  integral)  can  be  written  as  a  matrix 
operator.  In  general,  an  operator  can  be  represented  as  a  matrix  operator  with  respect 
to  a  set  of  basis  functions  [64],  For  this  convolution  integral,  the  matrix  operator 
representation  is  valid. 


Using  a  simple  one  dimensional  example  to  illustrate  the  method,  an  example  of 
converting  a  convolution  integral  to  a  convolution  sum  is  given.  Starting  with  a 
simple  form  (for  illustrative  purposes)  with  a  convolution  integral  [65] 


y{t)  =  |  h(t)x{t  -t)  dt 
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Now,  for  this  example  x(t)  is  a  general  time  series  and  h(t)  is  the  impulse  response  of 
a  linear  filter.  Equation  (4.9)  can  be  written  as  a  convolution  sum  for  discrete  data  as 

y{n)=  f^mh[n-k]  (4.10) 

k=-p 

where  n  is  the  time  index  and  k  is  the  delay  index.  To  show  that  this  can  be  written  as 

a  matrix  operator,  simply  write  Y=y[n]  n=T,2,3...  and  X=x[n]  n=l,2,3, . and  show 

that  for  windowed  (i.e.  finite  duration)  signals: 


Ml] 

+ 

1 

-s: 

i 

l 

•  h[\-(p) 

Ml] 

M2] 

= 

h[2-(-p) 

M2] 

M«]_ 

_h[n-(-p ) 

1 

i 

_M«] 

(4.11) 


By  carrying  out  the,  multiplication  it  is  seen  that  each  term  in  the  vector  y[n]  is 
obtained  from  the  matrix  multiplication  in  (4. 1 1)  (which  results  in  the  sum  4. 10)  for 
each  term.  So  in  matrix  notation  (4. 1 1)  can  be  written  as 

Y  =  hX  (4.12) 


By  an  analogous  derivation  for  the  2  dimensional  convolution  integral 


^(xi)' 

'£(*,, M)  •• 

K(x„x\) 

’W,(M) 

= 

KM)  *■ 

_K{Xj,x[)  •• 

(4.13) 
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Therefore,  the  integral  operator  can  be  written  as  the  matrix  operator  K(x,x’)  and  it 
follows  that  the  measurement  equation  is  in  the  correct  form  (4.2)  as 

=  £(*,*')  +  n^(x)  (4  14) 

This  completes  the  derivation  of  the  measurement  equation  for  our  wideband 
spreading  function.  Equation  (4. 14)  coupled  with  (4. 1)  enables  the  continuation  of 
the  innovations  approach  of  the  derivation  as  given  by  Kailath  [11], 


4.2.3  Kalman  Filter  Derivation  -  Innovations  Approach 

The  problem  is  posed  here  as 

Given  a  noisy  set  of  measurements  wm(i)  i=T,2,3...X  (X  being  the  current  state  X=[s 
t]  characterized  by  (4.13),  find  the  linear  minimum  (error)  variance  estimate  of  the 
state  characterized  by  (4.1).  That  is,  find  the  best  estimate  of  w0(X)  given  the 
measurement  data  up  to  time  (time  delay/time-scale  pair)  x  where  Wa(x)=[wm(l) 
wm(2) . wm(x ' )].  Here  Wm(x~)  is  the  set  of  observations  of  the  spreading  function. 

The  minimum  (error)  variance  criterion  is  given  by 


J{X\x)  =  E{w'0{X)w0{X)\Wm(x)} 


(4.15) 
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or  using  the  standard  conditional  mean  notation  and  noting  that  ~  denotes  estimate 

w0(X\  x)  =  E{w0(X)\1Vm(x)}  (4.16) 

which  states  that  the  estimate  of  the  true  spreading  function  at  the  “time”  given  data 
up  to  X  is  equal  to  the  expected  value  of  the  true  spreading  function  conditioned  on 
the  measurement  of  the  spreading  function.  Thus,  the  criterion  is  written  as 

J(X  |  x)  =  E{w'0(X  |  jc)  wD(X  |  x)  =  trace{P{X  \  x))  (4.17) 

Where  p  is  the  estimated  error  covariance. 

There  are  three  special  cases  of  the  state  estimation  problem: 

1 .  The  1-step  predicted  estimate. 

2.  The  filtered  estimate. 

3.  The  fixed  lag  smoothed  estimate. 

The  “filtered”  estimate  #  (jc  |  ;c)  is  of  interest  here  as  this  is  the  form  which  gives 
recursive  estimates  in  a  prediction/correction  form.  It  is  given  by 

wa(x  I  x)  =  E{w0(x)  I  Wm(x)}  J g) 

First,  the  dynamic  form  on  the  estimator  is  chosen  and  the  estimator  is  constrained  to 


be  linear. 
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In  the  dynamic  case,  the  linear  estimator  [51]  is 


Sw(*)‘  K,..WK,-'WE. -KmvWZ,  (4.19) 

where: 

(419a) 

Is  the  batch  minimum  variance  (MV)  spreading  function  estimate  and 

»  >  (4.19b) 

is  the  batch  measurement  and  Kmv  is  the  Kalman  gain. 

This  is  a  “batch”  minimum  variance  solution  to  the  state  estimation  problem,  since  all 
p-vector  data  {wm(l). . .  wm(N)}  are  processed  in  one  batch.  However,  for  the 
purposes  of  the  thesis,  a  recursive  solution  to  the  estimation  problem  of  the  form: 

A  A 

Wnew  =  Kid  +  KEnew  (4.20) 

is  needed. 

In  order  to  achieve  the  recursive  solution,  it  is  necessary  to  transform  the  covariance 
matrix  to  be  block  diagonal,  since 
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£K,(1K,(1)}  -  E{wm(i)wm  m 
E{wm(N)wJ(  1)}  E{wm(N)w^m 

Having  Ewm  be  block  diagonal  implies  that  all  of  the  off-diagonal  matrices 

RWy  ( i,j )  =  0  for  i^tj.  This  in  turn  implies  that  all  of  the  wm(x)  must  be  uncorrelated 

with  or  orthogonal  to  w'm(x). 


Therefore,  a  sequence  of  independent  p-vectors  {e(x)}  can  be  constructed  such  that 


E{e(x)e\k)}  =  0  forx  *  k 

The  sequence  e(x)  can  be  constructed  using  the  orthogonality  property 


(4.22) 


[wm(x)  -  E{wm(x)  \Wm(x-  \)}]lWm(x  - 1) 

Defining  the  innovation  or  new  information  as: 

e00  =  wm(*)-wo(x|x-l) 

with  the  orthogonality  property  that 

cov(wm(X),  e(x)]  =  0  forX  <x -l 


(4.23) 


(4.24) 


(4.25) 


Since  e(x)  is  a  time-uncorrelated  p-vector  sequence 
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Re(N)  = 


ReO)  - 

••  ReO) 

Re(!)  ■■ 

•  ReW 

(4.26) 


for  each  Re(i)  e  Rpxp. 


The  correlated  measurement  vector  can  be  transformed  into  an  uncorrelated 
innovation  vector  though  a  linear  transformation ,  say  L,  given  by, 


Wy  =  Le  L  e  Rpnxpn 


(4.27) 


where  L  is  a  non-singular  transformation  matrix  and  e:=[e'(l). . .  e'(N)] '. 


Multiplying  (4.26)  by  its  transpose  and  take  the  expected  value  to  obtain 


RWm(N)  =  LRe(N)L' 


(4.28) 


now,  inverting: 


(4.29) 


Similarly, 
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(4.30) 


Now,  substitute  these  results  into  (4. 19)  to  obtain 


LrW  =  K,»SK)Ky~\K)Ky  = 


(4.31) 


Simplifying  yields 

W^iN)  =  Kxe(K)K\K)e  (4.31a) 


Since  the  e(x)  are  time  uncorrelated,  Re(N)  is  block  diagonal.  From  the  orthogonality 
properties  of  the  innovations,  e,  it  can  be  shown  that  R*oe(0D  is  lower-block 
triangular,  that  is: 


Koe0D  = 


X  >  f 

K,e(x,x) 

■for- 

x  =  i 

0 

x<i 

(4.32) 


substituting  into  (4.31) 


wx(l,l) 

'Rv,e&  1) 

0 

X_,(  i) 

0 

wx(N,N)_ 

_Kxe(N,\) 

•••  R.xe(N,N) 

0 

r;1w 

e{N\ 

(4.33) 


Any  row  of  (4.32)  can  be  written  as 
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wx(x>x)  =  ZK„e(x>OK  \i)e(0 


i=\ 


If  the  last  term  of  the  sum  in  (4.33)  is  extracted,  from  (4. 18) 


(4.34) 


*,(■*,*)  =  £^(-*.0^0X0 + Kxe(x>x)K'(x)e(x) 


i=i 


W„, 


(4.35a) 


Or 


W„ew  =  ™x(*>*)  =  w(x  I  x  - 1)  +  K(x)e(x) 


(4.35b) 


Where  K(x)  =  RwAjc)(x, x)R;1  =  p(x \x-  1)C\x)r;\x)  is  the  Kalman  weight,  C  is  the 
measurement  matrix,  and  P(x  \  x  - 1)  is  the  updated  error  covariance  estimate. 


Recall  that 


™m(x  l*-0  =  K(x,x’)w0(x  |  x  - 1) 


(4.36) 


Using  the  measurement  equation  (4.14)  and  (4.24)  the  innovations  sequence  can  be 
decomposed  as 


e(x)  =  wm(x)  -  K(x,x')w0(x  |  x  - 1)  =  £(x,x')[w0(x)  -  wQ(x  \  x  - 1)]  +  W fn(x)  (4.37) 


e(x )  =  K(x,x')wa(x  |  x  - 1)  +  Wj-n(x) 


(4.38) 
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Consider  the  innovations  covariance  Re(x)  using  the  expression. 


Re  =K(x,x')P(x\ 


(4.39) 


and 

-  E{w0{x)e\x)}  =  £{w0(x)[/T(x,  x>c(x  |  x  - 1  )  +  W w(x)]'} 

ox  = 

or 


K,e  =  E{w0(x)w0(x  |  x  -  l)}£(x,  x') 


(4.40) 


Using  the  definition  of  and  substituting  for  wc 


K0e  =  E{[wo(x  |  X  - 1)  +  wa(x  I  x  - 1)]  w0(x  I  x  - 1)}£  (x,x')  (4.41) 


So  that 


K0e  =  £{[*«,(*  I  *  -  l)tf0(*  I  *  -  !)]}£(*,*')  +  E{[w0(x  I  x  -  \)wa(x  I X  -  !)]}£(*,*')] 


(4.41a) 


From  the  orthogonality  property  of  the  estimation  error  for  dynamic  variables,  that  is 


E{{Wm{X))w'0{x  |  jc  - 1)}  =  0  forX  <x-\  (442) 


The  first  term  in  (4.41)  is  zero  because  of  (4.42)  therefore  (4.41)  becomes 
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RWoe  =  P(x\x-l)K(x,x') 


(4.43) 


Therefore,  the  Kalman  gain  is 


K(x)  =  R„oe(X)R;l(x)  =  P(x\x- l)£(x,x’)R;l(x) 


(4.44) 


Before  the  corrected  state  estimate  can  be  calculated,  the  predicted  estimate  is  needed, 
that  is: 


WM  =  w0(x|x-l)  =  E{w0(x)  |  WJx  - 1)} 


(4.45) 


Employing  the  state  space  model  (4.1)  and  from  the  linearity  properties  of  the 
conditional  expectation, 

wx(x  |x-l)  =  E{A(x - 1  )wx(x - 1)}  (4.46) 


write  (4.46)  as 


^(x  |  x  - 1)  =  A(x  -  1)m^x  - 1 1  x  - 1) 


(4.47) 


or 


K,d  =  w0(x\x-i)  =  A(x-  l)*„(x  - 1 1  x  - 1)  (4.48) 


The  error  covariance  is  given  by 
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P(x  |  x - 1)  =  A(x - 1  )P(x -l\x- 1  )A'(x - 1)  (4.50) 

The  corrected  covariance  p(x  _  j  |  x  _  is  calculated  using  the  corrected  estimation 
error  and  (4.35a)  to  obtain 


w(x  I  x)  =  w0(x)  -  wQ(x  I  x)  =  w0(x)  -  wa(x I  AT  - 1)  -  K(x)e(x) 


or 


w(x  |  X)  =  w0(x  I  x  - 1)  -  K(x)e(x) 


(4.51a) 


(4.51b) 


Using  (4.51b)  the  required  covariance  is  calculated: 


P(x  |  x)  =  P(x  |  x  - 1)  -  K(x)E{e(x)w'0(x  I  *  - 1)}  -  E{w0(x  \  x  -  l)e'(x)}K'(x) 

+  K(x)Re(x)K'(x) 


(4.52) 


From  the  orthogonality  property  E{w0(x  \  x  -  l)e'(*)}  =  0  the  error  covariance  is 
now  given  by: 

P(x  |  x)  =  P(x  |  x  - 1)  -  K(x)K(x,x')P(x  |  jc  - 1)  -  P(x  |  x  -  1)£(jc,x')K'(x) 

+  K(x)^(jc)K’(*)  (4- 5 

Factoring  p(x  |  jc  - 1)  from  the  first  two  terms  and  using  the  expression  (4.44)  for 
K(x)  in  the  last  term 


P(x  \  x)  =  [/  -  K  (x)K(x,x')]P(x \x-l)-P(x\x- l)r(x,x')K\x) 

+  P(x  |  x  -  \)r(x,x')R;'Re(x)K'(x) 


(4.54) 
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Observing  that  the  last  two  terms  cancel  the  error  covariance  is  given  by 


?( x  |  x)  =  [I  -  K(x)K(x,x')]P(x  |  jc-1) 


(4.55) 


This  completes  the  derivation  of  the  Kalman  filter  algorithm  for  the  wideband 
spreading  function  recursive  updating. 


4.3  Analysis  of  Kalman  Filter  Derivation  Results 

In  this  section  the  major  results  of  the  previous  section  are  discussed.  First  the 
criterion,  the  mean  square  scattering  function  error  will  be  analyzed.  Second,  the 
dimensions  and  possible  uses  of  these  quantities  in  a  signal  processing  structure  will 
be  presented  and  an  explanation  of  how  the  quantities  are  interrelated  is  given.  Table 
4. 1  contains  the  results  of  the  Kalman  filter  derivation  for  recursive  scattering 
updating. 
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Table  4. 1  Summary  of  the  Kalman  Filter  for  Recursive  Scattering  Function 
Estimation  and  Updating. 


Prediction 

w0(x  |  x  - 1)  =  /t(x  -  i)wG(x  - 1)  +  q{x  - 1) 

(State  Prediction) 

P(x  |  x  - 1)  =  A(x  -  1)£(* -l|;c-  1M’(*  - 1) 

(Covariance  Prediction) 

Innovation 

e(x)  =  wm(x)  -wm(x\x-\)  =  wm(x)  -  K(x,  x')wjx  \  x  - 1] 

(Innovation) 

Re  =  K(x,  x')P(x  |  x  - 1)£  (x,  x')  +  RwgMx) 

(Innovation  Covariance) 

Gain 

k(x)  =  p(x\x-i)k;(x,x')r;1 

(Kalman  Gain/Weight) 

Correction 

w0(x  I  x)  =  wQ(x  1  x  - 1)  +  K(x)e(x) 

(State  Correction) 

P(x  |*)  =  [7-  K(x)K(x,x')]P(x  |  *  - 1) 

(Covariance  Correction) 

Initial  Conditions 

O 

O 

iST 

o 

o 

X 

The  rest  of  this  chapter  will  be  dedicated  to  the  analysis  of  the  quantities  outlined  in 
Table  4.1.1.  All  of  the  quantities  outlined,  however,  most  of  our  focus  will  be  on  the 
correction  terms,  the  innovations,  and  the  Kalman  gain. 
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4.3.1  Prediction  T  erms 

In  this  section  the  prediction  relationships  presented  in  Table  4. 1  are  discussed.  The 
state  prediction  is  simply  the  state  equation  (4. 1).  This  equation  gives  the  relationship 
between  the  values  of  the  predicted  spreading  function  at  the  current  (t,s)  conditioned 
on  (t,s)  at  the  last  time  t-1  (what,  in  the  derivation,  are  referred  to  as  x  and  x-1 
respectively).  The  predicted  spreading  function  at  these  two  “times”  are  related  by 
the  state  transition  matrix  A.  In  words,  the  new  estimate  of  the  spreading  function  is 
the  old  estimate  transformed  by  the  state  transition  matrix  A.  A  itself  is  derived 
directly  from  the  physics  of  the  problem  (See  Chapter  5).  The  state  transition  matrix 
in  the  WTD  related  the  new  scale  and  delay  to  the  old  is  based  on  the  motion  of 
objects  in  the  medium. 

The  second  item,  the  error  covariance  P ,  given  in  Table  4. 1  has  the  same  kind  of 
relationship  between  new  and  old  values  as  the  state  prediction  via  the  state  transition 
matrix.  In  this  case  the  new  value  of  the  error  covariance  is  the  old  value  transformed 
using  A  as  given  in  Table  4. 1 

The  assumptions  made  in  the  derivation  pertaining  to  the  prediction  parts  of  the 
algorithm  are  as  follows: 

1 .  A  is  a  known  matrix. 

2.  The  additive  process  noise  is  zero-mean  and  spectrally  white  with  correlation 


matrix  defined  by: 
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x  =  x' 
x±x' 


(4.56) 


4.3.2  Innovations 

The  innovations  approach  [55]  is  the  heart  of  the  derivation.  In  this  section,  the 
properties  of  the  innovations  or  “new  data”  [51]  and  the  innovations  relationship  to 
the  scattering  function  updating  algorithm  are  discussed.  Also,  the  innovations 
sequence  itself  ,e(x),  as  well  as  the  innovations  covariance,  Re(x)  will  be  discussed. 

The  first  property  used  in  the  derivation  of  the  Kalman  Scattering  Function  Algorithm 
is  that  the  Mxl  vector  e(x)  represents  the  new  information  in  the  observed  data  wy(x). 
The  innovations  process  has  the  following  properties  [66]: 

1 .  The  innovations  process  e(x),  associated  with  the  observed  data  Wy(x)  at  “time”  x 
is  orthogonal  to  all  past  observations  e(l)  e(2)...e(x-l)  given  by 


E{e(k)Wm(k)}  =  o  l<A<(x-l) 


(4.57) 


2.  The  innovations  process  e(x)  consists  of  a  sequence  of  vector  random  variables 
that  are  orthogonal  to  each  other: 


E{e{k)e'(k)}  =  0 


1  <  k  <  (x  -  1) 


(4.58) 
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3.  The  sequence  of  vector  random  variables  Wm(x).  and  the  sequence  of  vector 
random  variables  e(x)  have  a  one-to-one  correspondence.  This  means  that  one 
can  be  obtained  from  the  other  by  means  of  linear  stable  operators  without  the 
loss  of  information.  The  sequence  of  vector  random  variables  defining  the 
innovations  process  is  obtained  using  a  Gram-Schmidt  orthogonalization 
procedure  [67], 

As  seen  in  Table  4.1,  the  transition  “matrix”,  A,  and  the  measurement  “matrix”,  K  are 
actually  a  tensors.  These  matrices  can  be  transformed  to  a  more  recognizable  form 
size  by  realizing  what  these  tensor  notations  mean  and  by  using  a  lexiographic 
ordering  scheme  [68-70],  This  can  be  done  to  transform  the  tensor-matrix  form  to  a 
matrix  vector  form  without  the  loss  of  information.  However,  the  size  of  the  resulting 
transformation  matrix  can  be  problematic  for  implementation. 

In  the  next  section,  an  extended  Kalman  filter  which  tracks  the  peak  locations  within 
the  spreading  function  is  introduced  for  the  case  in  which  direct  measurements  of  the 
kinematic  quantities  are  not  available.  In  this  derivation  only  the  peak  locations 
within  the  spreading  function  are  tracked.  Time-scale  and  time-delay  are  estimated 
from  a  delay-scale  map.  However,  it  is  desired  to  track  state  variables  that  differ 
from  the  measurement  variables. 
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4.4  Extended  Kalman  Filter  Derivation 


In  this  section,  the  extended  Kalman  filter  [71]  to  track  a  single  highlight  moving  in  a 
circular  motion  is  derived.  The  state  variable  is  a  vector  quantity  containing  the 
velocity  vn  and  angle  0  as  shown  in  Figure  4.1. 


Scatterer 


O  Tx/Rx 

Figure  4. 1 :  Circular  Motion  Geometry 


It  is  assumed  that  the  highlight  has  an  initial  position  given  by  the  vector 


*0  = 


V0 

6, 


Which  are  uncertain  by  the  amounts  ctv  and  cr©.  The  true  highlight  “position”  moves 
randomly  about  its  unknown  nominal  location  with  standard  deviations  crv  and  gq. 
The  highlight  is  illuminated  with  a  pulse  and  the  time-scale  and  time-delay  of  the 
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return  are  measured  using  a  wideband  matched  filter.  This  measurement  has 
associated  with  it  error  due  noise  and  an  inherent  error  in  the  ability  of  the  transmitted 
pulse  to  estimate  the  true  time-delay  and  time-scale.  This  was  discussed  extensively 
in  Chapter  3  in  the  form  of  the  Cramer-Rao  lower  bound.  The  estimation  errors  on 
time-delay  and  time-scale  are  given  by  standard  deviations  ax  and  as.  The  derivation 
of  the  extended  Kalman  Filter  (EKF)  is  begun  by  first  introducing  the  state  equation 


xk  =  Axk,i  +  WkA 


(4.59) 


where  Wk_ f  is  the  state  noise  which  is  white,  Gaussian,  zero  mean,  and  has  a 
known  variance.  The  matrix  A  is  a  state  transmission  matrix  which  relates  past  to 
future  values  of  the  state  vector  which  is  given  as 


(4.60) 


Next  a  measurement  equation  is  needed  as  in  the  previous  sections,  however,  here, 
the  state  variables  cannot  be  observed  directly  as  the  measured  quantities  are 
functions  of  the  observed  quantities.  The  measurement  equation  is  given  by 


Z*  =Hk(xk)  +  Vk 


(4.61) 


where  zk  is  the  measurement  vector,  is  the  measurement  as  a  function 

of  the  state  variables  and  =[VS  Ve] '  the  measurement  noise,  which  is  white. 
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Gaussian,  zero  mean,  and  has  a  diagonal  covariance  R*  with  diagonal  elements 
defined  by  the  CRLB  derived  in  Chapter  3.  Our  measurement  equation  is  therefore 
given  in  terms  of  the  kinematic  variables,  and  observables  as 

c-v„(*) 

c  +  v„ 

—sin0(k) 
c 


0 


(4.62) 


Now,  to  start  the  iteration,  initial  estimates  of  the  state  and  the  error  covariance  are 
needed.  These  are  given  by 


1 

O 

1 _ 

cr2v„  0  1 

X 

© 

II 

1 

_ i 

n 

a? 

- 1 

<N 

b 

o 

_ 1 

(4.63) 


First,  this  initial  estimate  is  projected  ahead  to  the  first  measurement 


A 


(4.64) 


Now,  using  the  state  transition  matrix  derived  in  Section  5.3,  the  estimate  of  the 
variance  on  the  estimate  can  be  calculated.  This  is  given  as 


A(-) = a0p<>4 


(4.65) 


now,  the  update  of  the  state  variable  estimate  is  given  by 
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_  (4.66) 

which  when  4.62  and  4.64  are  substituted  becomes 


“a  ” 

V* 

An 

A. 

Ah 

+  Kk 


slk 

*k~ 


c  +  vk(~) 


c-v. 


2d 

—sin  0(~) 


(4.67) 


Now  the  Kalman  gain  Kk  is  calculated.  In  the  linear  case  it  is  a  function  of  Hk,  the 
measurement  matrix. 


Using  the  vector  Taylor  series  expansion 


H(x  +  Ax)  =  H(x)  + 


8H 

dx 


Ax  +  Higher  Order  Terms  (4.68) 


Having  made  the  assumption  of  negligible  acceleration,  the  higher  order  terms  are 
dropped.  This  is  only  possible  if  the  observation  interval  is  short  enough  to  ensure 
second  order  effects  are  negligible.  Dropping  the  higher  order  terms 


where 


AH  =  —&x 
dx 


(4.69) 


dH 

dx 


=  H 


(4.70) 
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and 


Are  the  elements  of  the  new  measurement  matrix  which  is  given  by 


H  = 


dHi 

3c, 

SH2 

3x, 


c'H, 

3c, 


dH,  dH, 

. .  •  •  •  — — — 

dx0  dx; 

: 

dx2 

•  • 

•  • 
t  • 

5H, 

•  «  •  •  •  •  

dxj 


(4.71) 


Which  in  our  case  is  a  2x2  matrix  for  the  time-scale  as  a  function  of  velocity  and  time 
delay  as  a  function  of  the  angle  0. 


The  iteration  can  now  be  completed  by  calculating  the  Kalman  gain  Kk 


Kk  =  Pk  (-)Hl  (xk  (~))[Hk  (xk  (~))Pk  (-)//[  (xk  (-))  +  Rk  I'1'4- 72> 


and  the  error  covariance  is  given  by 
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Pt(+)  =  [I-KkHt  (*(-))]/>*  (-)  (4.73) 

This  completes  the  iteration  of  the  extended  or  linearized  Kalman  filter.  This  is  the 
form  of  the  algorithm  used  in  the  remainder  of  the  thesis  to  implement  the  examples. 

4.5  Conclusions 

In  this  chapter,  a  Kalman  filter  for  updating  and  estimating  wideband  spreading 
functions  using  the  innovations  approach  has  been  calculated.  Through  the  use  of  the 
extended  or  linearized  Kalman  filter,  the  algorithm  that  tracks  the  kinematics  of  the 
problem  using  measurements  of  time-delay  and  time  scale  in  a  wideband  range-scale 
map  is  derived.  In  the  next  chapter,  the  extended  Kalman  filter  will  be  implemented 
for  objects  undergoing  various  types  of  motion.  The  ability  of  the  Kalman  filter  to 
maintain  a  track  on  the  peak  locations  within  the  spreading  function  of  the  object  will 
be  investigated  for  signals  with  varying  bandwidths  and  durations. 


Chapter  5 

Recursive  Kalman  Scattering  Function  Estimation  -  Implementation 


5.1  Introduction 


This  chapter  introduces,  via  two  examples,  the  implementation  of  the  recursive 
scattering  function  updating  algorithm.  The  first  example  will  be  of  two  rough  spheres 
attached  by  a  non-scattering  rod.  The  spheres  and  rod  rotate  as  a  single  object  in  the 
plane  which  contains  the  transmitter/receiver.  The  second  example  consists  of  a  set  of 
point  scatterers  linearly  arranged  and  moving  with  a  constant  velocity  and  with  only 
transverse  motion  with  respect  to  the  observer.  For  both  of  these  examples,  these  state 
and  measurement  Equations  (4.1)  and  (4.2)  will  be  calculated  explicitly,  and  the  effects 
of  transmitted  signal  parameters  on  KF  performance  will  be  presented  and  explored  in 
detail. 

5.2  Example  1:  Two  Rotating  Rough  Spheres 


For  this  example,  the  geometry  of  the  scatterers  is  as  shown  in  Figure  5.1. 


Figure  5.1:  Example  1  Geometry 


Observer 
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The  system  consists  of  two  rough  spheres  connected  by  a  non-reflecting  rod  which  is 
rotating  with  angular  frequency  oo.  The  system  will  be  illuminated  with  a  high  resolution 
VFM  pulse  with  various  durations  and  bandwidths  which  correspond  to  the  measurement 
error  variance  bounds  (CRLB)  given  in  Table  5.1. 


Table  5. 1 :  Transmitted  Signal  Parameters 


Signal 

Duration  (s) 

Bandwidth  (Hz) 

var[s  -  i] 

var[r  -  r] 

.25 

1000 

le-1 

le-1 

.5 

1000 

le-2 

le-1 

1 

1000 

le-3 

le-1 

1 

1 

le-1 

le-1 

VFM6 

1 

100 

le-1 

le-2 

VFM7 

1 

1000 

le-1 

le-3 

It  is  also  known  that  the  signal  to  noise  ratio  affects  the  CRLB  for  both  time-scale  and 
time  delay  as  discussed  in  Chapter  3.  However,  the  final  effect  on  the  Kalman  filter  is 
the  same  as  changing  the  transmitted  signal  parameters.  The  final  effect  being  that  the 
convergence  as  a  function  of  measurement  error  is  guided  by  the  error  covariance 
regardless  of  the  effects.  However,  it  should  be  noted  that  the  signal  to  noise  ratio’s 
effect  on  the  CRLB  is  multiplicative  as  discussed  in  Chapter  3,  and  therefore  is  an 
important  factor.  As  most  realizable  systems  are  reverberation  or  peak  transmitted 
power  limited,  it  is  the  signal  parameters  over  which  the  system  designer  has  the  best 
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control,  and  therefore  present  the  more  interesting  (and  useful)  results  of  the  Kalman 
filter  response  as  a  function  of  transmitted  signal  parameters. 

The  pulse  has  sufficient  time-bandwidth  product  to  resolve  the  rough  spheres  in  time  and 
frequency.  By  taking  the  auto-wavelet  transform  of  the  signal,  the  auto-ambiguity 
function  of  the  SFE  pulse  is  produced.  This  is  shown  in  Figure  5.3. 


Figure  5.3:  Auto-Ambiguity  Function  of  a  Is  1000Hz  VFM 

5.3  Derivation  of  the  State  Transition  Matrix 

In  this  section,  the  equation  of  state  for  the  system  presented  in  section  5.2  is  derived. 
By  considering  a  few  assumptions  about  the  system.  First,  neglecting  the  effects  of 
acceleration.  In  this,  it  is  assumed  that  the  system  is  rotating  slowly  enough  to  neglect 
angular  acceleration,  but  fast  enough  to  produce  a  measurable  scaling  of  the  return. 
Second,  a  reference  point  is  chosen  such  that  when  the  rod  connecting  the  spheres  is 
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perpendicular  to  the  line  of  sight  to  the  sensor,  the  relative  delay  is  zero.  The  scale  at  this 
point  is  a  maximum  due  to  the  fact  that  the  component  of  velocity  of  the  spheres  is 
parallel  to  the  line  of  sight,  making  it  a  maximum,  thus  producing  the  maximum  scale. 
When  the  rod  connecting  the  spheres  is  parallel  to  the  line  of  sight,  the  component  of 
tangential  velocity  parallel  to  the  line  of  sight  is  zero  -  producing  zero  scale.  Here  the 
delay  (+  and  -)  is  a  maximum.  Figure  5.4  shows  the  details. 


Figure  5.4  Scale  and  Delay  as  a  Function  of  Position  or  Time 


Now  knowing  each  scale  and  delay  as  a  function  of  time,  the  state  transition  matrix  is 
derived. 
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In  the  x-y  plane,  it  is  known  [72]  that  the  coordinate  transformation,  in  this  case  a 
rotation,  is  given  by  the  following: 


X 

cos# 

-sin# 

X 

_y_ 

new 

sin# 

cos  <9 

y_ 

old 

(5.1) 

where  9=0o+cot. 


The  trick  is  to  derive  a  similar  expression  for  the  new  scale  and  delay.  This  is  done  by 
writing  down  the  following  equation: 


-S3] 


old 


(5.2) 


Where  G  is  a  2x2  matrix.  The  following  are  also  used: 


1  + 


s—  =  — 7 ;  v  =  =-  cos(0o  +  cot) 
1--  2 


d 


(5.3a) 


2  d 

^w  =  -0‘+-sin(6'0+«O) 

c  2 


!  +  ■ 


Sold  ~ 


C  . 


l-V- 


■  v  =  -cos(0o) 


roid  =  ~(r +7sin(^0)) 
c  2 


(5.3b) 


(5.3c) 

(5.3d) 
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Using  (5.3a-d)  in  5.2  and  solving  for  G 


1  +  - 


( od  cos  (ait) 


2(c-<yycos(^0)) 

0 


1  +  - 


d  sin(<5tf) 


2(r  +  -sin(^0)) 


(5.4) 


-told 


Equation  (5.4)  is  the  form  needed  for  the  state  transition  matrix.  In  order  to  show  that 
this  is  the  correct  state  transition  matrix,  the  actual  and  noiseless  state  predicted  scales 
and  delays  are  plotted  in  Figure  5.4 


Figure  5.5:  Expected  and  Predicted  Values  of  Scale  and  Delay  vs.  Time 


Figure  5.5  shows  that  the  state  transition  matrix  produces  the  desired  results,  producing 
the  new  scale  and  delay  at  time  t,  given  the  measurement  at  time  t-1  -  meaning  the 


previous  measurement.  Now,  for  the  purposes  of  the  simulation  the  state  equation  is 
used: 
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™T\T-\(X)  -  A(tJ  -  l)w'j-_1(x) 


(5.5) 


This  gives  the  value  of  the  spreading  function  at  time  T  conditioned  on  the  previous 
value  at  T-l.  The  state  noise  is  considered  to  be  zero  mean,  white,  Gaussian  noise  with  a 
known  covariance  Rw. 


5.4  Rough  Rotating  Spheres  -  Data  Creation  Algorithm 

In  this  section,  the  algorithm  used  to  produce  the  data  necessary  for  the  examples  in 
this  chapter  is  discussed.  It  is  desired  to  track  time-scale  and  time-delay,  and  therefore 
two  data  sets  are  produced.  First,  the  time-scale  and  time-delay  for  the  model  portion  of 
the  algorithm  is  created.  Second,  the  measurements  of  the  actual  system  (i.e.  the  data) 
are  calculated.  Here,  the  noisy  measurements  are  produced  that  are  statistically  correct 
having  variances  of  white  Gaussian  noise  added  which  correspond  to  the  predicted  CRLB 
for  the  transmitted  signal  used.  Not  only  are  the  measurements  noise,  but  inaccuracies  in 
the  model  are  introduced  via  the  measurements.  The  idea  here  is  that  the  model  does  not 
accurately  represent  the  actual  system.  For  example,  in  the  examples  to  follow,  the 
distance,  d,  between  the  spheres  is  larger  for  the  simulated  measurement  data  that  for  the 
model.  In  this  way,  the  ability  of  the  Kalman  filter  to  track  measured  data  is  tested. 

Also,  the  point  is  found  at  which  inaccuracies  in  the  measurements  (as  a  function  of 
transmitted  signal  duration  and  bandwidth)  inhibit  the  ability  of  the  Kalman  filter  to 
estimate  the  actual  time-delay  and  time-scale  as  a  function  of  time.  Figures  5.6  and  5.7 


Figure  5.7:  Velocity  as  a  function  of  time  for  Measurement  and  Model 


Changing  only  the  distance  between  centers,  significant  differences  in  the  model  and 
measurement  quantities  can  be  produced  for  both  time-delay  and  time-scale  (theta  and 
velocity). 
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By  using  the  above  algorithm,  various  aspects  of  tracking,  detection,  and  signal 
design  can  be  related.  This  is  accomplished  by  observing  that  the  TW  improves 
processing  gain  [4]  (lowering  measurement  inaccuracies).  Work  in  Chapter  3  showed 
that  longer  durations  and  higher  bandwidths  lower  the  CRLB  for  estimates  that  can  be 
made  on  time-delay  and  time-scale.  This,  in  turn,  as  illustrated  in  the  following  section, 
affects  the  ability  of  the  Kalman  filter  to  track. 

5.5  Rotating  Spheres  -  Results 

The  results  of  the  processing  for  the  rotating  spheres  example  which  was  discussed 
above  are  now  given.  First,  the  results  for  time-delay  measurements  which  correspond  to 
the  angle  theta  (as  seen  in  Figure  5.8-5. 10)  are  given. 

It  is  desired  to  show  the  effect  of  changing  the  mean-squared  bandwidth  of  the 
transmitted  signal  on  the  Kalman  filter’s  ability  to  track.  Beginning  with  Figure  5.8 
which  contains  plots  of  three  quantities  vs.  time.  These  are  the  model  theta,  the  actual 
(measured)  theta,  and  the  updated  or  filtered  theta,  the  quantity  of  interest. 
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Figure  5.10:  Theta  Vs.  Time:  Transmitted  Signal  VTM  1  (Duration  Is,  BW  1  Hz) 

As  seen  from  the  above  series  of  plots,  the  bandwidth  of  the  signal  determines  the 
convergence.  Figure  5.8  shows  that  the  filtered  estimates  reflect  what  is  going  on  in  the 
actual  kinematics.  However,  as  the  bandwidth  is  decreased,  the  filter  follows  the  model 
completely  and  no  longer  tracks  the  actual  kinematics.  This  is  entirely  due  to  the 
effective  bandwidth  of  the  VFM  signal,  as  all  other  parameters  have  been  kept  constant. 
The  series  5. 5-5.7  is  interpreted  as  follows.  First  the  higher  bandwidth  VFM  has  enough 
resolution  to  accurately  track  the  spheres  as  they  rotate  as  seen  in  Figure  5.8.  Second, 
Figure  5.8  shows  that  a  lOx  reduction  in  bandwidth  causes  the  filter  to  diverge  (i.e.  ignore 
the  measurements).  Lastly  when  the  bandwidth  is  reduced  essentially  to  a  tone,  the  noise 
of  the  measurements  is  very  large  and  at  first  fluctuates  wildly,  until  the  gain  is  so  small 
that  the  measurements  are  ignored  completely. 
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5.6  Linear  Highlight  Motion  -  Results 

In  this  section,  the  results  of  Kalman  filtering  algorithm  operating  on  a  different 
kinematic  system  are  presented.  In  this  system,  a  highlight  in  linear  motion.  For  the 
simplicity  of  illustration,  the  results  for  the  single  highlight  are  given. 

5.6. 1  Linear  Motion  -  State  Transition  Matrix  Derivation 

In  this  section,  the  state  transition  matrix  for  a  single  highlight,  moving  with  constant 
velocity  v0  in  a  straight  line  is  derived.  The  geometry  for  this  scenario  is  given  in  Figure 
5.8. 


Figure  5.11.  Linear  Highlight  Motion  Geometry 
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The  state  transition  matrix  is  derived  as  before  in  Section  5.3.  However,  in  this  instance, 
the  normal  velocity  is  constant  and  therefore  An=l.  The  quantity  to  track  which  changes 
with  time  is  the  relative  delay.  This  is  given  by  the  standard  geometric/trigonometric 
relationships.  The  result  for  the  state  transition  matrix  is 


0 

v0  (t  + 1)  sin  <90 


(5.6) 


where  to  is  the  one  way  delay  time  from  the  transmitter/receiver  to  the  highlight  at  xq  and 
0o  is  the  angle  the  trajectory  the  highlight  makes  with  a  perpendicular  to  the  line  of  sight. 


This  example  differs  from  the  previous  example  in  that  the  extended  Kalman  filter  for 
the  delay  is  not  necessary.  Here  delay  appears  as  both  a  measured  and  modelled 
parameter.  This  changes  only  one  quantity  in  the  Kalman  filter  derivation,  which  is  the 
element  of  H22.  Here,  H22  =1  since  the  measured  quantity  is  simply  the  quantity  to  track. 
For  this  example,  all  of  the  parameters  used  in  the  first  example  are  also  used  in  this 
example.  The  signals  used  are  the  same,  with  changes  in  duration  and  in  bandwidth. 

The  algorithm  for  production  of  data  is  the  same,  however  the  parameter  changed  to 
produce  model  vs.  measured  data  is  the  angle  0o.  This  again  causes  the  modeled  velocity 
(normal)  and  delay  as  a  function  of  time  to  differ  from  the  measured. 
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5.6.2  Linear  Highlight  Motion,  Velocity  Tracking 

This  section  presents  the  results  for  the  velocity  tracking  portion  of  the  simulation. 
Three  results  are  presented.  First,  given  in  Figure  5.12  a  plot  containing  three  curves. 
These  are:  in  green,  the  filtered  output,  red  is  the  actual  trajectory,  blue  is  the  modeled 
normal  velocity.  In  this  case  the  measurement  accuracy  is  such  that  the  Kalman  filter 
tracks  immediately  to  the  data. 


5.1 

5 

4.9 

45 

4.7 

Velocity  4  a 
(m/s) 

45 

4.4 

45 

42 

4.1 

0  10  20  30  -40  SO  60  TO  80  90  100 

Time  (s) 

Figure  5.12:  Velocity  vs.  Time  Transmitted  Signal:  VFM  1 

The  next  plot.  Figure  5.13,  shows  how  the  filter  begins  to  lose  the  track  as  the 
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measurement  error  is  increased. 
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13:  Velocity  vs.  Time:  Transmitted  Signal:  VFM2 


plot  shows  how  the  filter  rejects  the  new  data  and  tracks  to  the  model. 
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5.6.3  Linear  Highlight  Motion:  Delay  (Range)  Tracking 

In  this  section  the  results  of  the  delay  portion  of  the  simulation  are  given.  Again,  the 
three  plots  given  correspond  to  the  signals  given  in  Table  5.2.  It  is  seen  that  as  the 
bandwidth  is  increased,  for  a  fixed  SNR,  the  ability  of  the  Kalman  filter  to  track  the  delay 
can  be  directly  related  to  the  transmitted  signal’s  bandwidth.  In  Figure  5. 14  it  is  seen  that 
the  Kalman  filter  tracks  to  the  received  data.  This  is  due  to  the  fact  that  the  measurement 
error  is  very  small,  and  therefore,  the  filter  “believes”  the  data. 
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Figure  5.15:  Range  vs.  Time,  Transmitted  Signal:  VFM4 
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Figure  5.16:  Range  vs.  Time,  Transmitted  Signal  VFM5 

Figure  5.16  shows  that  the  Kalman  filter  tracks  to  the  data  even  though  the  measurement 
errors  are  beginning  to  get  high.  The  filter  finally  breaks  down  and  begins  to  track  to  the 
model  in  Figure  5.17. 
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Figure  5. 17:  Range  vs.  Time,  Transmitted  Signal  VFM5 

5.7  Performance  Monitoring 

In  this  section,  a  method  of  determining  the  performance  of  the  Kalman  filter  is 
introduced.  This  method  includes  the  monitoring  of  the  spectrum  of  the  innovations 
sequence.  The  innovations  sequence  was  discussed  in  detail  in  Chapter  4.  It  is  known 
from  this  discussion,  that  for  the  optimal  estimator  [51]  the  innovations  sequence  is 
white.  In  this  section,  the  results  of  calculating  the  innovations  sequence  for  the  sphere 
moving  in  circular  motion  are  presented.  For  this  example,  plots  of  the  Kalman  filter 
output  for  the  velocity  and  for  the  angle  theta  are  given.  In  the  example  chosen,  the 
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From  the  above  Figure  5. 19,  it  is  seen  that  the  “whiteness”  of  the  spectrum  of  the 
innovations  sequence  is  a  good  indicator  as  to  the  operation  of  the  Kalman  filter.  One 
can  imagine  a  recursive  model  updating  algorithm  being  incorporated  into  the  state 
equation  portion  of  the  tracker  in  order  to  optimize  the  Kalman  filter  in  terms  of 
broadening  the  spectrum  of  the  innovations  sequence  and  thus  optimizing  the  Kalman 
filter. 


5.8  Conclusions 

The  above  observations  lead  to  the  conclusion  that  the  transmitted  signal’s 
parameters,  i.e.  duration,  bandwidth,  and  even  shape  (as  this  too  affects  the  CRLB  as 
shown  in  Chapter  3)  can  affect  the  ability  to  track  the  highlights  in  time-scale  and  time- 
delay.  It  has  been  shown  that  the  Kalman  tracking  algorithm  is  easily  reconfigured  for 
different  Kinematic  scenarios.  The  main  point  of  this  work  is  to  show  that  the  issues  of 
signal  design,  detection,  and  tracking  can  be  intimately  related.  The  above  analysis 
shows  that  the  theoretical  calculations  of  the  Cramer  -Rao  lower  bound  on  the  variance 
of  the  estimates  that  can  be  made  have  a  direct  impact  on  the  Kalman  filters  ability  to 
converge  to  the  actual  scenario  and  to  balance  the  model  and  measurement  portions  of  the 
tracking  algorithm.  The  above  examples  have  shown  that  a  model  based 
detection/tracking  algorithm  can  be  designed  for  which  there  exist  observables  (i.e., 
innovations  and  Kalman  Gain)  which  provide  a  direct  measure  of  the  performance  of  the 


tracker. 


CHAPTER  6 


Summary,  Conclusions  and  Further  Research 

6.1  Summary  of  Results 

The  primary  focus  of  this  thesis  has  been  to  derive,  implement  and  investigate  a  recursive 
algorithm,  based  on  a  model  based  Kalman  filtering  approach,  to  estimate  and  update  the 
location  of  peaks  in  the  wideband  spreading  function  which  correspond  to  the  kinematic 
properties  of  an  object.  In  this,  primarily  interest  has  been  in  the  physical  variables  of 
position  and  velocity  with  the  actual  observables  being  time-scale  and  time-delay  in  a 
range-scale  map  produced  by  a  wideband  matched  filter.  Several  original  contributions 
were  made  in  the  completion  of  this  task. 

The  primary  contributions  made  in  this  thesis  are: 

•  Derivation  of  the  multidimensional  Kalman  filter  for  total  spreading  function 
updating  and  estimation. 

•  Reduction  of  the  dimensionality  of  the  above  derivation  to  handle  updating 
and  estimation  of  the  peak  locations  in  a  range-scale  map. 

•  Relating  the  state  variables  to  the  measured  quantities  in  the  extended  Kalman 
filter  derivation. 

•  Derivation  of  the  Cramer  -  Rao  lower  bound  for  time-scale  and  time-delay  for 

the  wideband  signal  model  and  wideband  detector. 

•  Implementation  of  the  Kalman  filter  for  spreading  function  updating  and 
executing  illustrative  examples. 
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•  Provide  a  theoretical  and  experimental  (through  simulation)  common  thread 
which  begins  with  the  transmitted  signal  design,  ambiguity  properties  of  the 
transmitted  signal  and  the  estimation  properties  of  the  signal  including  the 
effects  of  noise  (CRLB)  and  connects  these  properties  to  the  eventual 
convergence  issues  of  the  extended  or  linearized  Kalman  filter. 

The  thesis  began  by  introducing  the  necessary  background  for  development  of  the 
Kalman  filtering  algorithm.  This  included  a  cursory  review  of  wavelet  transforms, 
Kalman  filtering,  and  applications  of  the  output  of  our  filter  to  model  based  signal 
processors  in  the  form  of  the  wavelet  transform  estimator-correlator. 

Then,  using  the  slowly  fluctuating  point  target  signal  model,  the  Cramer-Rao  lower 

* 

bound  on  the  estimation  errors  for  time-scale  and  time-delay  for  wideband  transmitted 
signals  was  derived.  This  contribution  is  important  as  this  relates  transmitted  signal 
properties  such  as  alpha  moment,  envelope  shape  in  the  time-frequency  plane,  bandwidth, 
and  duration,  to  variances  on  errors  associated  with  making  estimates  on  time-delay  and 
time-scale  using  a  wideband  matched  filter  detector. 

Having  the  relationships  discussed  in  the  previous  paragraph,  the  Kalman  filter  which 
tracks  the  kinematic  properties  of  the  system  based  on  the  measurements  (estimates) 
made  on  time-delay  and  time-scale  using  the  wideband  probing  signal  can  be  discussed. 
Within  the  derivation  of  the  Kalman  filter  a  measurement  equation  becomes  an  important 
quantity.  In  this  case  a  matrix  form  of  the  wideband  matched  filter  is  the  measurement. 
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Associated  with  the  measurement  equation,  is  a  quantity  referred  to  as  the  measurement 
error,  this  measurement  error  is  assumed  to  be  zero  mean,  Gaussian,  and  white  with  a 
known  variance,  namely,  the  Cramer-Rao  lower  bound  is  used.  In  this  way,  the 
transmitted  signal  parameters  are  related  to  the  errors  associated  with  the  measurements, 
these  errors  can  then  be  related  to  the  ability  to  track  the  kinematic  parameters  of  the 
physical  system  using  the  Kalman  filter.  It  should  also  be  mentioned  that  the  parameters 
of  the  detection  process,  such  as  signal  to  noise  ratio  also  appear  in  the  Cramer-Rao  lower 
bound. 

Following  the  derivation  of  the  total  wideband  spreading  function  predictor-corrector 
version  of  the  Kalman  filter,  the  problem  is  reduced  form  a  tensor-matrix  problem  to  a 
matrix-vector  problem  by  using  a  lexiographic  ordering  technique.  In  this  way,  an 
extended  Kalman  filter  that  tracks  only  the  highlights  of  interest  (peaks  in  the  range-scale 
map  which  are  local  phenomena  for  which  the  CRLB  is  valid).  In  this  way  state 
equations  which  predict  the  kinematics  of  multiple  highlights  can  be  derived,  thus, 
reducing  the  dimensionality  of  the  problem.  In  this  way,  a  model  updated  state  equation 
is  used  with  a  wideband  matched  filter  measurement  equation  to  predict  kinematic 
properties  based  on  the  past  measurements. 

Once  the  extended  Kalman  filtering  algorithm  was  derived,  a  computer 
implementation/simulation  was  made.  Using  the  simulation,  two  data  sets  were 
produced.  The  first  data  set  corresponded  to  the  model.  The  second  set  corresponded  to 
the  physical  system.  These  were  intentionally  made  to  be  different  to  test  the  ability  of 
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the  Kalman  filter  to  track  the  physical  system  even  though  the  model  did  not  accurately 
reflect  the  physical  processes.  In  other  words,  the  ability  of  the  Kalman  filter  to  balance 
the  accuracy  of  the  measured  data  to  the  expected  values  predicted  by  the  model  (state 
equation)  was  tested.  In  this  way,  it  could  be  shown  that  the  signals  with  wider 
bandwidth  produces  lower  CRLB’s  on  time-delay,  and  thus  allowed  us  to  have  better 
estimates  on  time-delay.  This  in  turn  allowed  the  Kalman  filter  to  better  track  the 
physical  position  of  the  scatterer.  Also,  it  was  shown  that  longer  duration  signals 
produced  lower  CRLB’s  on  time-scale.  This  in  turn  showed  that  the  accuracy  of 
estimates  on  time-scale  could  be  better  “believed”  by  the  Kalman  filter.  This  in  turn  was 
used  to  show  that  the  Kalman  filter  could  better  track  the  true  velocity  of  the  scatterer. 
Interestingly,  plotting  curves  of  the  CRLB  for  delay  and  scale  showed  that  signal  to  noise 
ratio  and  duration  or  bandwidth  could  be  “traded  off’  in  the  sense  that  if  a  lower  SNR 
exists,  a  higher  bandwidth  or  longer  duration  signal  can  be  used  to  overcome  the 
shortcoming.  Figure  6. 1  illustrated  the  topics  covered  and  their  interrelation.  Using  the 
above  framework,  each  of  the  blocks  in  Figure  6. 1  illustrates  a  starting  point  for  further 
research.  This  is  the  topic  of  the  next  section. 
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Figure  6. 1 :  Outline  of  Topics  Covered  and  Their  Interrelationship 

6.2  Recommendations  For  Further  Research 

In  this  thesis  several  questions  which  are  stand  alone  topics  for  further  research  have 
arisen.  In  this  section,  avenues  for  pursuing  these  topics  are  discussed. 

In  Chapter  3,  CRLB’s  for  time-delay  and  time-scale  were  derived.  Examples  for  the 
uncoupled  estimated  by  making  assumptions  on  the  shape  of  the  signal,  and  by  using  a 
real  signal  were  produced.  One  important  topic  for  further  research  is  to  do  an  exhaustive 
study  of  the  full  CRLB  for  any  signal.  This  will  involve  the  more  difficult  to  interpret 
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cross  term  J12.  As  this  is  an  important  signal  design  issue,  an  extensive  parametric  study 
should  be  completed  to  show  how  the  transmitted  signal  parameters  affect  coupled 
estimates. 

Next,  the  Kalman  filter  should  be  implemented  for  a  larger  number  of  highlights  for  a 
more  realistic  object  should  be  investigated.  This  should  include  an  exhaustive  study  of 
the  tuning  parameters  for  tracking  the  highlight  structure  of  a  complex  object.  Included 
in  this  study  should  be  a  study  of  the  robustness  of  the  approach,  as  well  as  a  study  of  the 
effects  of  errors  in  the  initial  estimates  should  be  investigated. 

Stability  issues  should  also  be  investigated  along  with  the  issue  of  efficient 
implementation  for  a  larger  number  of  highlights.  Bierman  [73]  has  extensively 
investigated  the  problem  of  efficient  implementations  for  Kalman  filtering  routines. 

Lastly,  the  cases  for  non-white  and/or  non-Gaussian  measurement  and  state  noises  should 
be  investigated.  In  this  case,  the  noise  in  the  measurements  should  include  reverberation 
terms.  This  in  itself  is  an  important  and  self  contained  research  topic,  namely,  deriving 
the  CRLB  for  time-delay  and  time-scale  under  reverberation  limited  conditions. 

In  the  larger  picture  it  would  be  desirable  to  further,  in  a  fundamental  sense,  the 
multidimensional  issues  raised  in  Chapter  4.  Remembering  that  in  Chapter  4,  it  was 
shown  that  the  state  equation  and  measurement  equations  contain  four  dimensional 
matrix  multiplications.  Mathematics,  specifically  tensor  analysis  coupled  with  operator 
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theory  should  be  investigated  deeply  to  gain  better  understanding  of  the  relationships 
between  these  tensor/matrix  operations  and  operator  theory.  In  the  case  of  the  state 
equation,  it  was  shown  that  each  sub-matrix  operating  on  a  state  variable  was  a  state 
transition  of  that  variable.  The  tensor  relationship  in  the  measurement  equation  is  far  less 
clear.  Usually,  to  dispense  with  the  higher  dimensionality,  assumptions  such  as  wide- 
sense-stationary  -uncorrelated  scattering  are  invoked.  Better  mathematical 
understanding  of  these  mechanisms  may  allow  the  relaxation  of  some  of  these 
assumptions  and  provide  a  more  general  solution  to  the  wideband  problems  at  hand. 

Further,  work  completed  here  should  be  used  to  highlight  the  importance  algorithms  for 
multi-feature/highlight  object  tracking.  Future  research  should  investigate  the  use  of 
learning  algorithms  for  application  to  model  based  signal  processing.  This  would  enable 
systems  to  learn  from  the  environment  and  change  their  own  operating  parameters  to 
better  interact  with  the  environment.  These  would  be  truly  intelligent  systems.  In  the 
future,  autonomous  systems  should  be  developed  to  apply  learning  algorithms  to  update 
statistical  models.  Applications  would  include  image  tracking,  dynamic,  real-time, 
medical  imaging  displays,  high  resolution  radar  and  sonar,  and  intelligent  vision  systems 
for  industrial  robots.  The  ideas  begun  in  this  thesis,  if  taken  to  their  ultimate 
development  would  enable  robotic  systems  to  create  a  “mental  image”  of  an  object  and 
track  its  motion  based  on  the  robot’s  “knowledge”  of  kinematics  (i.e.,  range  of  possible 
motion).  Using  the  kinematics  in  combination  with  the  ability  to  follow  object  features 
through  rotations  and  translations  (i.e.,  various  aspects  over  time)  provides  a  robust 
tracking  ability  for  future  detection/estimation  systems. 
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APPENDIX  A 


Derivation  of  the  Jn  Term  of  the  Fisher  Information  Matrix 


In  this  section  a  derivation  of  the  Ju  term  of  the  information  matrix  is  provided  in 
order  to  illustrate  the  methods  used  in  the  derivation  of  all  of  the  components. 


Beginning  with  the  equation  for  the  J 1 1  term  from  Chapter  3 


2c  Re 


dL  dL* 
dr  dr 


+  L 


d2L 

dr2 


(A!) 


and  the  detector 


L(s ,  t)  =  J  r(t)4sf*  ( s(t  -  r))dt 


(A.2) 


also  needed  is  the  echo  model  for  r(t)  as  given  by 


r(t)  =  b^[U(x)f(t)]  +  n(t) 


(A.3) 
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where  Et  is  the  transmitted  energy  and  b  is  a  random  variable,  white,  Gaussian 
distributed,  and  zero  mean.  Also,  n(t)  is  a  zero-mean, white,  Gaussian  randon  variable. 


Here 


[V(x)f(t)\  =  ^f(s(l  ~  T» 

(A.4) 


Substituting  A(3)into  (A.2),  inserting  the  results  into  (A.  1)  and  taking  the  prescribed 
derivatives 


2  c 


2alE,\\[U(x)f(t)W{x)m}’ 


OT 


-§-[U(x)f(u)]dtdu 

or 


[2c jj  N0S(t  -  )]*  ~^[U(x)f(u)]dt  du 

or  or 


(A.5) 


+ 


2c 


2cr2bEt  J f  [U(x)f(t)][U(x)f(u)]\U(x)f(  0]*  ±-[U(x)f(u)]dt  du 


dr 


2c 


N0  JJ  S(t  -  u)[U(x)f(u)]*  ±-[U(x)f(u)]dt  du 


dr 
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Now,  to  reduce  the  first  term  in  (A.  8),  we  use  Parsevals  Relation  to  change  to  the 
frequency  domain 

]s[sf(s(t-T))e-JMdt 

—CO 

Using  the  substitutions. 


2  =  J  J  4sf(s{t-r))e~Ja*dt  j- 

-oo  -oo  lUT 


x  =  s(t-  t) 


(A.  10) 
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Integrate  Equation  (A.  10)  to  obtain 


(A.  11) 


Which  is, 


(A.  12) 


Now,  turning  our  attention  to  the  second  term,  it  is  seen  that  following  the  same 
procedure,  that  is  using  Parsevals  relation  and  making  the  same  substitution  using  (A.  10) 

Jn  -  CO2  -CO2  (A.  13) 

Which  is  the  same  as  for  the  narrow  band  case,  this  tells  us  that  the  CRLB  for  time  scale 
is  inversely  proportional  to  the  mean  squared  bandwidth. 
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